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Plasma instabilities can be encountered in many branches of physics. This work 
focuses on relativistic plasmas with applications in theoretical astrophysics and 
particle accelerator physics. Even though these fields seem to be unrelated the un- 
derlying plasma physics processes arc often very similar. Two plasma instabilities 
- the beam-beam instability and the coherent synchrotron radiation instability - 
are analyzed. The former severely limits the achievable luminosity in storage rings 
and is related to the two-stream instability which has been proposed as a candidate 
for the radiation mechanism of radio pulsars. The main emphasis is on coherent 
synchrotron radiation which can lead to prohibitive energy losses in bunch com- 
pressors. Coherent synchrotron radiation also makes up the intense emission of 
radio waves by pulsars. Simple models based on the linearized Vlasov equation 
and relativistic magnetohydro dynamics which allow to compute detailed spectra 
of the emitted radiation are developed. 
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PREFACE 



Plasma physics effects are ubiquitous. Most of space is made of plasma and many 
devices on earth (including household appliances) generate plasmas, e.g. fluores- 
cent bulbs, klystrons in microwave ovens, electron beams in CRTs, particle accel- 
erators or electron microscopes etc. There are two recurring main questions which 
arise in plasma physics: How do plasmas evolve in time and are there equilibria 
which are stable under the influence of small perturbations? How much energy 
is lost due to radiation? Even though the mentioned applications do not seem 
to have much in common, plasma instabilities and radiation are often caused by 
the same few well-known (and some not so well-known) processes. This opens up 
possibilities for testing astrophysical processes in the laboratory. 

After covering a few basics, the beam-beam instability (a rather unpleasant 
instability encountered in storage rings which severely limits the achievable lu- 
minosity) is reviewed. In some aspects this instability resembles the two-stream 
instability which is currently considered to be responsible for the radio emission 
of spinning neutron stars (radio pulsars). Chapter [5] deals with the Coherent Syn- 
chrotron Radiation instability as an alternative to the two-stream instability in 
radio pulsars. According to the preceding paragraph it may not come as a surprise 
that this instability also plays a crucial role in particle accelerators. In chapter [7] 
a simplified approach is presented which is based on magnetohydrodynamics and 
the results of a computer simulation thereof can be found in chapter [61 

It is hoped that one day a gifted experimenter will exploit these similarities 
and come up with particle accelerator experiments which might greatly benefit the 
astrophysics community. 
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Chapter 1 
Plasma Physics 

1.1 Statistical Mechanics of Plasmas 

In classical mechanics the trajectory of a particle x(t) is completely determined if 
the forces acting on the particle are known and if the position xq and the velocity 
Xq (or alternatively the momentum po) at an arbitrary initial time are known. 
One needs two initial conditions because the equations of motion are second order 
differential equations excluding some peculiar special cases. For a huge number of 
particles, e.g. gas molecules in a steel cylinder, it is obviously not very practical 
to compute the trajectories of all those particles, and it is not even necessary. 
Since for all practical purposes the gas molecules are indistinguishable even for a 
"classical" observer one can ask instead how many particles f{'x,p,t)d^xd^p one 
could encounter at time t, position x and momentum p inside a phase space interval 
of size d^xd^p assuming the normalization 

d^x j d^pf{x,p,t) = N (1.1) 

where is the total number of particles. Of course any other normalization is 
equally good. If all the forces acting on this collection of particles are known it 
should be possible to find an operator acting on the distribution function / which 
describes its time evolution. Such an "operator" does exist and the resulting equa- 
tion is known as the Vlasov equation. Note that a statistical treatment based 
on a smooth distribution function Q| eliminates certain features known as discrete 
particle effects which can have a rather big impact on the distribution function. 



^In simple models of shock formation the distribution function may tend to a 
discontinuous limit. 
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Recovering these effects is usually rather difficult but nevertheless important. In 
plasma physics for instance, continuous charge distributions like a rotating ring 
with completely uniform charge density do not radiate, i.e. even synchrotron radi- 
ation (both coherent and incoherent) is due to the discreteness of electric charge. 

1.2 Vlasov Equation 

Consider the following "microscopic distribution function" 

TV 

F(x, P,t) = Y, 5[x - x,(t)]5[p - p,(t)] (1.2) 

which contains all the information about each individual particle. It satisfies the 
following equation of motion 

t {I + V, ■ A + Ip, . A} .[X - - pM = (1,3) 

j=l 

(Klimontovich equation). In a plasma the average force exerted on a particle by all 
other particles (which acts like an external force) is bigger than the force exerted 
by the nearest neighbors. With 

{F) = f 

(p.) = P • (1-4) 
averaging the Klimontovich equation gives the Vlasov equation 

The last equation in (11. 4p is justified if binary correlations are assumed small. 
Thus, all discrete particle effects like binary collisions have been removed. They 
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can be recovered by writing Eq. (11. 4p as a sum of an average (external) part and an 
internal part due to nearest neighbor interactions. In section 11.81 this will lead us 
to the Fokker-Planck equation. For astrophysical plasmas the densities are usually 
so low that the plasma can be assumed to be coUisionless. 

1.3 Solving the Vlasov Equation 

For particles interacting electromagnetically one replaces in Eq. fll.51) with the 
Lorentz force 

^p = g(E + vxB) (1.6) 
The fields are related to the sources by the Maxwell equations 



V-B = (1.7) 

V-E = 47rp (1.8) 

VxE = -|b (1.9) 

d 

V X B = 47rj + — E (1.10) 



where 



p = qj rfW(x,p,t) (1.11) 

3 = qj rfV/(x,p,t) (1.12) 



Furthermore, velocities and momenta are related by 

p/m 



1.13) 



a/1 + p^/m^ 

where m is the rest mass of the particles. This system of equations is referred to as 
the relativistic Vlasov-Maxwell system of equations. In the purely electrostatic case 
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the Maxwell equations simplify to the Poisson equation and the resulting system 
is known as the Vlasov-Poisson system. The system of equations is nonlinear in 
/ and is therefore hard to solve analytically. Eq. fll.l2p can usually be linearized 
either in the non-relativistic or in the ultrarelativistic case. Finding an equilibrium 
distribution, i.e. a distribution with df/dt = 0, simplifies the system and this 
task is often doable under more or less realistic assumptions. It can be shown 
that any distribution function / which only depends on the constants of motion 
represents an equilibrium. Typically, the constants of motion are the Hamiltonian 
and the canonical angular momentum. However, rewriting / as a function of x 
and p still requires solving a nonlinear equation. Once an equilibrium is found one 
can ask whether this equilibrium is stable to small amplitude perturbations. This 
question can be answered by linearizing the Vlasov equation about the equilibrium 
distribution. A variety of instabilities has been analyzed this way. 

1.4 Solving the Linearized Vlasov Equation Using the Method 
of Characteristics 

As was pointed out earlier the nonlinear nature of the Vlasov equation makes it 
hard to find analytical solutions for it, but looking for equilibrium solutions may 
be a successful endeavour in many simplified situations. Once an equilibrium has 
been found it is possible to linearize the Vlasov equation about the equilibrium. 
This provides valuable information about the stability properties of the found 
equilibrium. Writing the distribution function and the fields as a sum of the 



equilibrium part and a perturbation 

f = f + Sf 
E = E° + 

B = B° + 5B (1.14) 
and neglecting second order terms one obtains the linearized Vlasov equation 

{ I + - ■ - » + ^ X =°) ■ ^} W ^ -e (^E + V X B) . ^/» (1.15) 

The linearized Vlasov equation can be rewritten by following a particle on an equi- 
librium orbit (x', p') which passes through (x, p) at time t' = t. The equilibrium 
orbits have to satisfy the equations of motion 

Ax'(0 = v'(0 (1.16) 

^p'it') = e {E"(x'(t')) + v'(t') X B"(x'(t'))} (1.17) 

Thus, 

^Sfi^'it'), p'(t'),t') = -e (5E(x', t') +v'x <5B(x', t')) ■ (x', p') (1.18) 
Integrating the last equation an expression for Sf can now be obtained easily. 

1.5 Mathematical Properties of the Vlasov- Meixwell Sys- 
tem 

Since it is hard to develop an intuition for how solutions of the Vlasov-Maxwell 
system look like for a particular set of initial data it is desirable to find theorems 
regarding general properties of the system. Will the solution remain smooth at all 
times for smooth initial data, i.e. will the momentum space carrier of the distribu- 
tion function remain compact in finite time? Will the solution remain symmetric 
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for symmetric initial data? For the three dimensional relativistic Vlasov-Maxwell 
equations the answer to the first question is still unknown, but lots of incremental 
progress has been made. A very good review article on this subject is |2]. Glassey 
and Schaeffer [3] have shown that the answer is "yes" in two spatial and three 
momentum dimensions. The answer is important for deciding whether shocks can 
form spontaneously or not. The answer to the second question is "yes" for spheri- 
cally symmetric initial data [2]. 

1.6 Magnetohydrodynamics 

The distribution function / contains all the meaningful information one could 
possibly ask for in a statistical treatment and it is determined by either the Vlasov 
equation in the absence of discrete particle effects or by the Boltzmann or Fokker- 
Planck equation in the presence of discrete particle effects. These equations are 
usually difficult to solve and solving them numerically is usually not an option 
either - at least not in three dimensions where one would have to deal with seven- 
dimensional PDEs. Even though numerical solutions are being obtained by Ellison 
and collaborators [4J for lower dimensional problems (one spatial dimension, two 
dimensions in phase space) this is usually not a viable method. 

The underlying problem is that the distribution function still contains a lot of 
information and one may wonder whether the problem simplifies if one is content 
with less information. For many practical purposes it may be enough to compute 
certain macroscopic quantities like 

the number density 
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the mean velocity 



u(x,t) 



= n 



1 



v/(x,p,t)rf^p , 



(1.20) 



the mean momentum 



p(x,t) 



= n 



I 



P/(x,p,t)rf^p , 



(1.21) 



and the three dimensional stress tensor 



P = m / (v - u) (g) (v - u)/(x, p, t)d^p 



(1.22) 



i.e. moments of the distribution function. The definition of the stress tensor 
depends on the previous definition of the mean velocity. 

More macroscopic quantities can be obtained by multiplying the integrand by 
an arbitrary power of momentum and / or velocity components. Starting from 
the Vlasov equation one can find equations determining those moments easily. 
Integrating Eq. (11.51) over all momenta gives the continuity equation. 



Multiplying Eq. (11.51) by v and integrating over all momenta results in the 
"Euler equation" 



These are exactly the same equations one encounters in fluid dynamics if the 
higher order moments, i.e. the stress tensor, is neglected. Therefore, this approach 
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is known as magnetohydrodynamics (MHD). An arbitrary number of equations can 
be found this way. There is one serious problem, though. The moment expansion 
does not close in the absence of collisions, i.e. each equation will couple to the 
next higher moment. In fluid dynamics the equation of state can be used to close 
the system, but in plasma physics such an equation is generally unavailable - at 
least in the absence of collisions. However, under certain conditions it may be 
possible to guess an equation of state or it may be possible to neglect the higher 
order moments, e.g. for a cold plasma in the absence of a pressure gradient and 
heat flux. Magnetohydrodynamics will usually give good results if the frequency 
which is characteristic for the evolution of the distribution is much smaller than 
the plasma frequency and the cyclotron frequency. 

1.7 Some Plasma and Fluid Instabilities 
1.7.1 Negative Mass Instability 

A longitudinal bunching can occur in a beam executing circular motion if the ef- 
fective mass of the particles is negative, i.e. if an increase in energy leads to a 
decrease in angular velocity {d<p/dE < where <p = |0|). The energy at which 
the sign of d(j)/dE changes is called transition energy. In a weak focusing machine 
(e.g. charged particles moving perpendicular to an external magnetic field with- 
out gradient) this condition is always satisfied. Assume that due to an arbitrary 
initial perturbation the charge density is higher at a certain point on the circle. 
The electrostatic potential tries to repel particles away from the center of higher 
density. Particles in front of the region of higher density gain energy, but their 
angular velocity decreases. Similarly, particles behind the region of higher den- 
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sity lose energy and increase their angular velocity. Thus, neighboring particles 
are attracted to the region of higher density. The negative mass instability can 
be compensated by a sufficiently large energy spread. In the limit of zero energy 
spread the dispersion relation is \5\ 

Nre rj 1 



27rro 7'^ (Ati;)^ 

where 



;i.25) 



6 dp 



and 



^ (1.27) 

N is the number of electrons, vq is the radius of the orbit, is the classical electron 
radius, and the azimuthal mode number is denoted by m. Classical references are 

mmu- 

1.7.2 Rayleigh- Taylor Instability 

The Rayleigh Taylor instability is a fluid instability which can develop if a less 
dense fluid with density pi propagates in a denser fluid with density p2- Clumps 
of gas observed in supernova remnants are often due to this instability. 

1.7.3 Kelvin-Helmholtz Instability 

The Kelvin-Helmholtz instability is a non-relativistic fluid instability which can 
form at the interface of two flows with different velocities Ui and U2. The ubiquitous 
water waves caused by wind blowing over the surface of a pond are a typical 
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example. The Kelvin-Helmholtz instability may also be important to understand 
the observed patterns of astrophysical jets surrounded by the interstellar medium. 

The complex frequency of a perturbation in a system unstable to both the 
Rayleigh- Taylor instability and the Kelvin-Helmholtz instability is given by [8] 



where k is the wavenumber of the perturbation and g is the gravitational acceler- 
ation. 

1.7.4 Diocotron Instability 

The Diocotron instability is ubiquitous in the circular motion of a low density non- 



devices. This electrostatic instability resembles the Kelvin-Helmholtz instability 
in the sense that it forms in the presence of shear. For equilibrium configurations 
with d/dz = a. sufficient condition for stability is that the number density n(r) 
is a monotonically decreasing function, i.e. the maximum number density occurs 
at r = 0. 

1.7.5 Cyclotron Maser Instability 

A relativistic beam moving along a guiding magnetic field may be subject to the 
cyclotron maser instability. Unlike the classical Diocotron or the negative mass in- 
stability it is a transverse electromagnetic instability which is capable of producing 
coherent electromagnetic waves. The instability is driven by an inverted popula- 
tion in the transverse (i.e. parallel to the magnetic field) momentum distribution, 
i.e. the momentum distribution is sharp and its mean is non-zero. The cyclotron 





neutral plasma with d<p/dr ^ and can be found in common microwave generating 
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maser instability is exploited in microwave generating devices. 
1.7.6 Two-stream Instability 

The two-stream instability is an electrostatic instability which occurs for a plasma 
consisting of two (or more streams) of (not necessarily the same species of) particles 
with different velocities. Its development requires a region in phase space with 
df/dp>0 and the momentum of the particles satisfying the latter inequality have 
to be large compared with the momenta of the remaining particles. This instability 
is of importance for the understanding of problems associated with the solar wind. 
In particle accelerators secondary emission of electrons from the beam pipe may 
provide a background of electrons which can trigger a two-stream instability. 

1.8 The Fokker-Planck equation 

The averaging employed in section 11.21 was rather crude and removed all discrete 
particle effects. In this section it is attempted to recover the effect of statistical 
processes like radiation damping and quantum excitation occurring at random 
times ti. Instead of Eqs. fll.4p one now uses 



The probability for the occurrence of a perturbation in momentum space with 
Axi and Api is given by the probability densities Px{Axi) and Pp{Api), respec- 
tively. Px and Pp are taken to be normalized to unity and symmetric in their 




(1.29) 
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arguments. Without the former condition the distribution function would not re- 
main normahzed. After a timestep At the phase space element AxAp changes by 
the factor 

as can be seen by expanding the evolution of AxAp to first order in At using 
Eq. fll.29p and (11.301) . Assuming the number of particles is conserved one obtains 

f(x + gAt,p + hAt,t + At)AxAp 1 + ( ^ + ^ ) At = 

L \ox op J 

/OO POO 
/ d{Ax)d{Ap)P.,:{Ax)Pp{Ap)f{x- Ax,p- Ap,t) (1.32) 
-oo J —oo 

Expanding the distribution function inside the integral to second order in Ax and 
Ap allows one to evaluate the integral. Making use of the properties of and Pp 
mentioned above 



dt dx dp \dx dp J 2 ^ dx"^ 2 ^ dp"^ 



where the coefficients and Hp are related to the second order moments of Px and 
Pp, respectively. Eq. (11.33P is known as the Fokker-Planck equation. It describes 
the evolution of a plasma under the additional infiuences of radiation damping 
due to incoherent synchrotron radiation (in this case the parenthesis in Eq. (11.331) 
differs from unity) and quantum excitation due to the statistical nature of the 
radiation process (the plasma emits "discrete" photons). The former tends to 
increase the phase space density whereas the latter tends to decrease it. Setting 
the parenthesis to unity and neglecting the excitation coefficients and Sp the 
Vlasov equation is recovered. 
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1.9 A Simple Solution of the Fokker-Planck Equation in 
Beam Physics 

Since radiation damping and quantum excitation counteract each other equihbria 
may exist, i.e. distribution functions with df/dt = 0. In action angle variables 
(cf. §2.41) one can find the following 0-independent equilibrium 

f-^^e-i (1.34) 

Beams whose equilibrium distribution is given by Eq. fll.34p are called Gaussian 
beams. They can be encountered in electron-positron rings with significant syn- 
chrotron radiation. 



Chapter 2 

Physics of Particle Accelerators 

2.1 Applications and Limitations 

Particle accelerators have become an invaluable tool for high energy physics exper- 
iments. Due to the increasing complexity of such machines particle accelerators 
themselves have become the subject of detailed theoretical studies. Current ma- 
chines can reach center of mass energies of up to 2 TeV (Tevatron / Fermilab, 
as of 2001) and luminosities of up to 6350 pb^^ per year (CESR / Cornell, as 
of 2000). Such parameters give rise to all kinds of instabilities. One can divide 
accelerators into circular and linear machines. The best known linear accelerator 
(LINAC) is probably the Cathode Ray Tube which accelerates electrons emitted 
from a filament by a large electrical potential difference between the filament and 
a plate with a hole in it. Due to the large potential differences needed for a high 
energy beam such accelerators become technically unfeasible beyond lOMV. In a 
Wideroe LINAC alternating current instead of direct current is used to acceler- 
ate the particles. Charged particles travel through an array of conducting tubes 
with alternating polarity. Negatively charged particles in a gap between two tubes 
which leave a tube at negative potential are attracted by the next tube at positive 
potential. When the AC source reverses the polarity the particles are inside a tube 
and are shielded from the fields exerted by neighbouring tubes. To account for the 
increasing velocity of the particles the tubes have to increase in length. Linear ac- 
celerators hke the proposed Linear CoUider tend to be rather long if high energies 
are desired. In a synchrotron the particles execute circular motion and can pass 
the accelerating structure multiple times before reaching the desired energy and 
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being injected into a storage ring, for instance. The increasing energy of particles 
passing through the accelerating structure can be taken into account by adjusting 
the frequency of the voltage applied to the structure. As the name suggests the 
purpose of a storage ring is to store the accelerated particles. Usually there is a 
rotating and a counter-rotating beam consisting of particles and antiparticles, re- 
spectively Q At the interaction point (there may be multiple ones) the beams cross 
each other and colliding particles may annihilate and produce new particles which 
can be detected by a huge detector surrounding the interacting point. However, 
the probability for such an event to happen is small and huge amounts of energy 
would be wasted if the remaining particles were just dumped. The idea of the stor- 
age ring is to keep the particles circulating in the ring (possibly for many hours) 
until they finally collide. There are two main problems with circular machines, 
though. Accelerated charges (in this case the acceleration stems from forcing the 
particles onto a circular orbit) emit synchrotron radiation. The energy loss due 
to synchrotron radiation can make the operation of such machines prohibitive at 
high energies. At ultra-relativistic energies the radiated power is given by 

n = i-a.^ (2.1) 

where is Sand's radiation constant for electrons 

47r T 

C = — — " = 8.8575 • lO^^mGeV-^ , (2.2) 
3 [meC^y^ 

Tc is the classical radius of the electron, R is the radius of the ring and E is the 
energy of the electron. Thus, for high energies very large radii are needed which 
makes accelerators expensive to build. Even if the energy loss is not a concern the 



^The never-completed SSC was supposed to collide protons onto protons (in- 
stead of antiprotons) 
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highest achievable energy is hmited by the magnet technology. Currently, the high- 
est fields are provided by superconducting magnets, but superconductivity breaks 
down at sufficiently high magnetic field strengths. The current record for a contin- 
uous field is 45.1 Tesla measured at the National High Magnetic Field Laboratory 
at Florida State University [H] . The simplest conceivable circular accelerators con- 
sists only of dipole magnets which bend the beam and two plates with a potential 
difference which accelerate the beam. Such an accelerator is known as a weak- 
focusing machine. Its drawback is the beam size which increases with increasing 
radius. Bigger and bigger machines were built until the apertures of the magnets 
became prohibitively big and expensive to produce. 



2.2 Strong Focusing 

The beam size was drastically reduced once strong-focusing machines were in- 
vented. These machines contain quadrupole magnets with alternating gradients 
(in addition to the dipole magnets). Therefore, the beam size and the force due to 
the quadrupole magnets depends on the position s in the ring which is a number 
between zero and the circumference 2ttR of the ring. The focusing force in the 
horizontal and vertical plane, respectively, is 

x" = -K^{s)x y" = ~Ky{s)y (2.3) 

with 

1 1 



KJs) = — + 



i?2 B^yR dx 



where x and y are the horizontal and vertical displacement from the design orbit, 
respectively. Derivatives with respect to s are denoted by a prime. Let us focus 
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Figure 2.1: Geometry of the strong focusing machine 

on the equation for the horizontal motion. A solution which satisfies the initial 
conditions x(0) = xq, x'(0) = x'q, w{0) = wq, w'{0) = w'q and ipi^) = i'o is given 

by 

X = [{wq sin'?/'o + 'Wq'^ cosipo) Xq — Wqx'q sin-^/'o] w{s) cos tpi^s) + 

[— (wq cos-j/^o — '>^o'^ sin-j/'o) xq + x'qWq cos^^o] w{s) sm'ip{s) (2.5) 

where the width w{s) of the beam is determined by the envelope equation 

w{s)" + K{s)w{s) + — 1-^ = . (2.6) 

Furthermore, 

ij'{s) = w{s)-^ . (2.7) 

Individual particles oscillate ("betatron oscillations") about the design trajectory 
i> times, but all particle orbits are contained in the "envelope" whose width is given 
by w{s). The w""^ term in Eq. (12.61) acts like a "centrifugal barrier" and gives the 
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beam envelope a non-zero width, v is called the (machine) tune which is defined 



as 



u= — (h (3{s)K{s)ds (2i 



271 

where 

Pis) = w'is) (2.9) 
is called the betatron function and f denotes the integration from zero to 27ri?. 

2.3 Weak Focusing 

In the case of weak focusing, i.e. K{s) = K, the above equations simplify dramat- 
ically. One obtains 

K = I3-'^ iy = l , (2.10) 

i.e. even in the absence of quadrupole magnets the beam executes one betatron 
oscillation per revolution. 

2.4 Emittance 

The forces exerted by dipoles, quadrupole and higher order magnets which may be 
needed to correct for certain "optical errors" are conservative, i.e. the phase space 
density occupied by the particles in a beam is constant. The phase space area in 
(x, x') space is vre where e is called the emittance. Rewriting Eq. fl2.5p as 

x{s) = v^e/3(s) cos/xo (2-11) 

the emittance is determined by 

e = -foxl + 2aoXoXQ + Pox'^^ (2.12) 
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where 



7 



(2.13) 



(2.14) 



(2.15) 



Similarly, for x' 





(sin fiQ — a cos /io) 



(2.16) 



Instead of using x and x' to describe the motion of a particle trajectory it is very 
often advantageous to use the so-called "action angle variables" ip J = e/2 
because the latter is a constant of motion if the system is conservative. Note 
that particles in an accelerator are subject to many non-conservative forces like 
synchrotron radiation damping and acceleration. 

2.5 Beam-Beam Interaction 

All the limitations mentioned in section 12.11 are well known, but there is a myr- 
iad of less obvious issues that arise from the collective behavior of the particles 
which interact electromagnetically among themselves. The most severe limit on 
the achievable luminosity is due to the beam-beam interaction. If two oppositely 
charged beams which are slightly off-axis collide head-on the rotating beam is de- 
flected by the electromagnetic field of the counter-rotating beam and vice versa. 
The beam-beam force is highly non-linear. Its presence causes a tune shift ^ at 
the interaction point. It is customary to estimate ^ using the linear part of the 
beam-beam force. ^ is proportional to the number density of the beam and there- 
fore ^ is frequently used to parametrize the strength of a beam. In e~^e~ colliders 
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the observed beam-beam limit is in the range 0.02 < < 0.1 |12l [13]. A further 
increase in ^ increases the vertical emittance and leads to particle loss. Attempts 
to cancel the beam-beam force have failed so far. In the DCI experiment at SAL, 
Orsay, France, pairs of electron and positron beams were made to collide, i.e. both 
beams had zero net charge [H]. The result was disappointing. No significant im- 
provement of the beam-beam limit was observed. This outcome was explained by 
Derbenev in terms of a collective instability of the four-beam system [15]. There- 
fore, it is reasonable to assume that a collective instability is responsible for the 
beam-beam instability in the two-beam system as well. Indeed, collective oscilla- 
tions are seen in computer simulations [151 [IS [U] • The linearized Vlasov equation 
has been used to study the stability of colliding beams dHl [IT]. In [T7] Chao 
and Ruth analyzed the stability properties of beams that are confined to motion in 
the vertical direction ("flat beams"). They perturbed a "water-bag" equilibrium 
which has a uniform density within an ellipse in the phase-space {y,y'). However, 
electron and positron beams tend to have a "Gaussian" distribution (cf. section 
II. 9p . Very roughly speaking there are more particles in the inner region of the 
ellipse than in the outer region. In chapter H] the stability properties of an electron 
beam colliding head-on with a positron beam are investigated. They beams are 
assumed to be fiat and have a "Gaussian" equilibrium distribution. Both angular 
and radial modes are considered. Radial modes are modes which change the size 
of the ellipse. It is found that the radial modes have a profound influence on the 
stability of the system. 



21 

2.6 Coherent Synchrotron Radiation 

The beam-beam interaction is ubiquitous in storage rings when beams colhde, but 
it is by no means the only significant instabihty. More recently an instability due 
to coherent synchrotron radiation (CSR) is being thoroughly investigated which 
has been identified as a potential problem for the design of the proposed linear col- 
lider. Since the particles in the beams have only one chance to collide one would 
like to decrease their emittances as much as possible in order to achieve high lumi- 
nosity. Such low emittance beams can only be produced in damping rings where 
the emittance is reduced by emitting synchrotron radiation. The linear collider 
needs very short bunches to operate, but beams in a damping rings are subject to 
other instabilities if their bunch length is too short. The solution is to reduce the 
bunch length in a device known as a bunch compressor before the beam is injected 
into the linear collider. Bunch compressors consist of an accelerating section and 
an arc section. In the arc the beam emits synchrotron radiation whose wavelength 
may be close to the bunch length. In this scenario the electromagnetic waves can 
modulate the beam in such a way that the bunches are equidistant. The radiation 
from individual bunches can now interfere constructively and the incoherent radi- 
ation becomes coherent. For coherent radiation the total radiated power scales as 
N"^ instead of where N is the number of particles. The beam would lose all its 
energy almost instantly. Therefore, it is important to know under which operating 
conditions one can avoid this effect. CSR has been observed already in a couple 
of accelerator labs [TSl [12] . In chapter [S] a simple model of a collisionless, rela- 
tivistic, finite-strength, cylindrical layer of charged particles is presented which is 
capable of emitting coherent radiation. The particles interact with their retarded 
electromagnetic self-fields in a way that allows them to clump together. Including 
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the radial dynamics is difficult, and a small energy spread (which translates into 
a small non-zero thickness of the rotating layer) is one of the main requirements. 
It is shown that the betatron oscillations can lead to a significant decoherence 
which is responsible for the emission of a very characteristic spectrum. The sta- 
bility properties are analyzed by solving the linearized Vlasov-Maxwell system of 
equations. The treatment resembles work by Uhm, Davidson and Petillo [20] who 
examined the stability of a thin relativistic electron ring. However, their interest 
is in the negative mass instability and their approximations are not suitable for 
electromagnetic effects like CSR. A simpler model in which all particles were con- 
strained to move on a circle with fixed radius was presented in 1971 by Goldreich 
and Keeley [21j. In this model the particles initially move at constant speed, but 
they can gain or lose energy by interacting with the azimuthal component of the 
electric field. However, no mechanism for fixing the radial degree of freedom is 
provided. It is unclear whether (or under which conditions) the radial degree of 
freedom can be neglected. If, for instance, the circular motion is due to an external 
magnetic field without gradient ( "weak focusing" ) an increase in energy translates 
into an increase of the orbit radius whereas the velocity remains almost constant 
in the case of ultra-relativistic motion. This may not be very favorable for the 
development of a bunching instability. A more realistic model was investigated by 
Heifets and Stupakov. In f22l 123] they analyze the stability of electrons executing 
circular motion. The radius of the individual particle orbit is determined by the 
energy of that particle, i.e. the radial motion and the relative longitudinal motion 
are coupled such that the problem has effectively only one degree of freedom. It is 
not entirely obvious under which conditions such an approach is valid. The model 
was extended by Byrd [24j to include the effect of a conducting beam pipe which 
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can serve as a cut-off of the allowed wavelengths. A conducting beam pipe can 
severely attenuate the CSR instability. 



Chapter 3 

Physics of Rotating Neutron Stars 

3.1 Stellar Evolution 

A living star is supported against its own weight by the pressure it builds up as a 
result of heat generated in fusion reactions inside the star. A young star generates 
its heat from the conversion of hydrogen into helium by nuclear fusion. Once the 
supply of hydrogen is exhausted in the core the star starts to shrink increasing 
its temperature. This allows the star to burn the remaining hydrogen in its shell. 
One has to distinguish two cases. 

3.1.1 M < 8Mq 

In the red giant stage the shell expands leaving behind the core which continues to 
shrink until a white dwarf is formed. Fluid instabilities destroy the shell turning 
it into a nebula. This stage can be regarded as the end point of the evolution 
of a light star. The white dwarf continues to emit thermal radiation until it has 
completely cooled down. 

3.1.2 M > 8Mq 

Heavier stars become hotter during contraction triggering fusion reactions of heav- 
ier elements. Fusion stops once all material in the core has been converted into 
ironB Like in the previous case burning continues in the outer shell. Instead of a 
red giant a super red giant is formed with a radius bigger than 100 million kilome- 



^The element with the highest binding energy is ^^Ni and not ^^Fe. Cf. an 
article by Fewell [22] on why iron is more abundant than nickel. 
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ters. The core is supported by the degeneracy pressure of non-relativistic electrons 
and - as the star continues to contract - the electrons become relativistic and the 
increase in pressure slows down. Furthermore, at relativistic electron energies the 
protons can capture electrons which turns them into neutrons, thus reducing the 
degeneracy pressure of degenerate electrons. Photodissociation of iron leads to a 
polytropic index smaller than 4/3 rendering the core unstable to collapse [26j. The 
iron core implodes which generates a shock wave propagating outward. The shock 
wave comes to a stop before it can leave the super red giant. However, under cer- 
tain conditions a bubble can form between the core and the shock front. A small 
fraction of the binding energy of the star is used to eject all the material of the 
star except the core in a supernova explosion. Its mechanism is complicated, but 
it is believed to be caused by convection and neutrinos transporting energy. The 
remnant is called a neutron star because it is only supported by the degeneracy 
pressure of degenerate neutrons. Sometimes the conditions under which a super- 
nova explosion takes place are not satisfied or an insufficient amount of matter is 
released. In this case even the degenerate neutrons cannot prevent the star from 
collapsing even further and a black hole is formed. 

3.2 Properties of Rotating Neutron Stars 

The radius of a typical neutron star is in the order of i? ~ 10 km while its mass 
is in the order of Mq. In addition to a strong gravitational field on their surface 
they also posses a strong magnetic field which can be as strong as 10^ Tesla. 
The field can be described by a magnetic dipole to a good approximation (cf. 
section 13. 3p . Like on earth the magnetic field may be created by the dynamo 
effect. Since charged particles are necessary to create a magnetic field a neutron 
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star cannot consist entirely of neutrons. Indeed, it is believed that a neutron star 
contains a small fraction of electrons and protons in its core [27] . Neutron stars 
born with a large amount of angular momentum are capable of emitting intense 
electromagnetic radiation in the radio frequency range. Such radio pulsars have 
rotation periods ranging from Is down to 33ms for the Crab pulsar. In general 
the axis of rotation does not coincide with the alignment of the magnetic dipole 
moment. Therefore, the radiation sweeps out a cone about the axis of rotation. 
Every time the observer's line of sight coincides with the magnetic axis a pulse of 
intense electromagnetic radiation is observed with a period equal to the period of 
the rotation of the star T = 27t/Q (lighthouse model). 

The discovery of millisecond pulsars ruled out white dwarfs as possible candi- 
dates for radio pulsars as can be seen from the following simple argument. For 
a stable star the centrifugal force exerted on particles at the surface of the star 
cannot exceed the force due to gravity. Thus, 



White dwarfs are not sufficiently dense to satisfy the inequality above. 

Due to its large mass and angular momentum the star posses a huge amount of 
kinetic energy which powers the emission of the intense radiation. Therefore, as the 
pulsar continues to lose energy its angular velocity has to decrease ("spin-down"). 
The details of the process converting kinetic energy to electromagnetic radiation 
are not completely understood yet. However, an argument based on conservation 
of energy suffices to relate some fundamental parameters of a pulsar. Assuming 
the energy loss is due to magnetic dipole radiation 



n^R < GMR 



(3.1) 



dW 



p = Wr'^b', 



(3.2) 



dt 
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Figure 3.1: Radio pulsars are spinning neutron stars with strong magnetic fields. 
At the magnetic poles their radiation follows the magnetic field lines. Since the 
field lines and the axis of rotation are misaligned the radiation sweeps out a cone. 
The region in which the radiation is believed to be generated is shown in gray. A 
similar cone could be drawn on the other side of the star. 
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i.e. measuring the spin-down Q and equating the energy lost by a magnetic dipole 
to the loss of kinetic energy one can solve for the normal component Bq of the 
magnetic field at the magnetic pole if the angular velocity Q, the radius R and the 
moment of inertia / are given. The spin-down can be measured very precisely. 

3.3 Braking Index 

Eq. (I32D relates n to (l 

QocQ'' (3.3) 

where the so-called braking index is denoted by n. According to Eq. (13.21) the 
braking index is 3 if the emission is due to dipole radiation. Deviations from n = 3 
would suggest that the simple model leading to Eq. (13.21) is not completely accurate. 
Indeed, braking indices as low as = 1.4 have been measured. The determination 
of the braking index is very simple if Q is known. Differentiation of Eq. (13.31) and 
making use of the same equation again to eliminate the proportionality constant 
one obtains 

n = ^ 3.4 

3.4 Some Fundamental Parameters of the Crab Pulsar 

We start by estimating how many charged particles could be present in the magne- 
tosphere. In a model by Goldreich and Julian [28] the axis of rotation is assumed 
to coincide with the orientation of the magnetic dipole moment. Assuming the 
neutron star and its surrounding magnetosphere along the magnetic field lines are 
perfect conductors one obtains 

E + (rixr)xB = (3.5) 
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Figure 3.2: The observed spectrum of the Crab pulsar extends from the radio 
regime to frequencies up to lO^^Hz pJ.The straight hne in the radio regime is 
proportional to u'^/^ 

Thus, the magnetosphere must have the Goldreich- Julian charge density 

ncj = (47r)-V • E = ~ lQ^^cm-\B /lQ^^G){R/rf[T /Is]-^ (3.6) 

at radius r > R where T = 27r/|f2|. 

The energy loss can be determined from Eq. fl3.2p by measuring Q and Cl. 

P = An'^lf/T^ (3.7) 
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It is in the order of 10'^^ erg s^^ for the Crab pulsar where [29] 

M ~ lO^^kg R ~ lO^m / ~ lO^^kg 
T = 33ms t = 4.22 ■ lO^^^ b = 5.2 ■ lO^^G (3.8) 

The highest detected frequency of the Crab pulsar is in the order of lO^^Hz, 
but the spectrum in Fig. 13.21 starts to drop off significantly at lO^^Hz. 

Integrating the Goldreich- Julian charge density from the surface of the surface 
to the velocity of light cylinder (Fig. 13. 3p gives 



N = 10^^cm-%B/10^^G)[T/ls]'^R^ An ■ r^drr'^ 

Jr 

= 10^^cm-%B/10^^G)[T/ls]~^ -AnRHn^ (3.9) 
~ 10^^ . (3.10) 



3.5 Emission Mechanism 

Obviously, one would like to have a better understanding of how the rotational 
energy is converted into radiation and in particular how the radiation mechanism 
works. In this paragraph it is shown that incoherent synchrotron radiation cannot 
account for the observed brightness of the radio signal. The synchrotron radiation 
is partly reabsorbed by the inverse Compton effect (cf. section [3. lip . The ratio of 
the brightness temperature (cf. section I5.10p due to inverse Compton radiation to 
the brightness temperature due to synchrotron radiation is given by 

1 



l + |(T^a../10^'K)5(/,/lMHz) 



(3.11) 
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For Tmax < lO^^K and an upper cutoff frequency fc ~ lO^MHz in tlie radio regime 
this ratio is smaller than one, but for Tmax > lO^^K 

Lic/L, ~ (T^,,/10"K)i° (3.12) 

Therefore, brightness temperatures exceeding lO^^K are impossible to achieve with 
incoherent synchrotron radiation (P oc A^) and some sort of coherent radiation 
mechanism (P oc A^^) is required. The brightness temperature of the Crab pulsar 
is roughly lO^^K. 

3.6 Secondary Electron-Positron Plasma 

The strong magnetic field forces the electrons to move parallel to the field. Since 
magnetic fields cannot do any work a strong electric field with E-B ^ is necessary. 
Because Eq. ( 13. 5 p implies E-B = the accelerating field must be due to a deviation 
from the Goldreich- Julian charge density [31] • Several effects accomplishing this 
have been suggested, e.g. general relativistic effects [32] or the bending of the 
magnetic field lines |33j . 

Some photons emitted by the accelerated charges create secondary electron- 
positron pairs which screen the electric field except in compact regions called 
"gaps". The particles are accelerated in those gaps. The plasma consisting of 
secondary particles has a distribution which differs significantly from the distribu- 
tion of primary particles. Instabilities of the primary plasma lead to the coherent 
emission of radio waves whereas instabilities of the secondary plasma lead to a 
non-thermal emission in the high frequency regime from IR to 7-rays. This work 
focuses on instabilities found in the primary plasma (which may be induced by the 
interaction with the secondary plasma). 
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^ cQl . 

Figure 3.3: Geometry of the regions surrounding a neutron star. The closed mag- 
netosphere is followed by a gap in which strong electric fields accelerate charged 
particles. The star is surrounded by a co-rotating magnetosphere which cannot 
extend beyond the velocity-of-light cylinder, i.e. the radius at which the velocity 
of the particles would exceed the speed of light at angular velocity fl where fl is 
the angular velocity of the star (and the co- rotating magnetosphere). 
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3.7 Free Electron Maser Emission 

Several mechanism were proposed to explain the emission of electromagnetic waves 
by the primary plasma. This paragraph deals with the so-called free electron maser 
emission. It requires a strongly modulated electric field parallel to the magnetic 
field. How such an electric field could be generated in space is unknown. Rowe 
|34] found that such a set-up is capable of self-amplification if the distribution is 
inversely populated, i.e. there is a region in phase space where the particle density 
increases as the energy increases. The electric field accelerates the particles which 
therefore emit electromagnetic radiation. This radiation then modulates the beam 
until bunches of particles radiate in phase (Actual free electron lasers which are 
built in labs use a magnetic field from an undulator instead of an initially modu- 
lated electric field. The undulator causes the particles to move on a helical path 
which then emit synchrotron radiation.). Unfortunately, the growth rate as a func- 
tion of energy decreases too rapidly to explain the high brightness temperatures 
that are observed. 

3.8 Two-stream Instability 

Particles of the secondary plasma might interact with those from the primary 
plasma and (due to their very different distributions in phase space) trigger a two- 
stream instability. All conditions for the development of a two-stream instability 
are met. However, detailed calculations [35] show that the expected growth rates 
are too low. Again, the high Lorentz factors and the low density of the involved 
beams are the offending parameters. Another problem is posed by the inability 
of the waves generated by the two-stream instability to escape from the neutron 
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star, i.e. they have to be converted into different waves which can actually escape. 
This may involve some yet unknown non-linear effects. Despite these shortcomings 
the two-stream instability is considered to be the most promising candidate for an 
explanation of the observed brightness temperatures by many authors [36J. Two- 
stream instabilities due to electrons streaming against positrons in the secondary 
plasma were also considered by several authors, e.g [37]. Again the reader is 
referred to the review article by Usov [36j. 

3.9 Curvature Radiation 

Curvature radiation was the first emission process which was studied in the context 
of radio pulsars. The radiation is due to the synchrotron radiation emitted by 
accelerated charged particles where the acceleration originates from forcing the 
particles to move along an arc. Since the observed radio waves are polarized it is 
natural to attribute them to synchrotron radiation. The coherence was explained 
by a maser-like mechanism by many authors. Most approaches based on maser 
curvature emission ran into trouble because the conditions for a self-amplifying 
maser instability were not satisfied. On the other hand non maser-like mechanisms 
which started out with a bunched distribution were heavily criticised because it 
was unclear how the bunches could form and because of a lack of detailed models 
which took the velocity spread of the distribution into account. 

In chapter one such model is presented. It is assumed that the radio emission 
is due to coherent curvature radiation which is produced by small bunches of 
particles whose radiation interferes constructively. For this approach the linear 
stability properties of a cylindrical, coUisionless, relativistic layer made of charged 
particles whose axis of rotation is aligned with an external magnetic field are 
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analyzed using the linearized Vlasov- Maxwell system. The particles are allowed 
to interact with their own electromagnetic self-fields. The bunches are seeded by 
arbitrarily small initial perturbations which grow exponentially in time until the 
perturbations saturate. Knowledge of the saturation amplitude is a prerequisite for 
calculating the intensity of the electromagnetic radiation and it can be estimated 
considering the trapping of particles in the "potential well of the wave" . 

3.10 Beaming 

If an isotropic emitter moves at relativistic speed an observer at rest observes 
the radiation as if it was radiated into a narrow cone pointing into the forward 
direction. Its opening angle is approximately F"^ where F is the Lorentz factor of 
the moving source. This effectively increases the power measured by an observer 
who can sample only a small solid angle. 



where a prime denotes quantities measured by the observer. 

Furthermore, the emitted frequency undergoes a relativistic Doppler shift. For 
a source moving towards the observer [38] 




3.11 Inverse Compton Radiation 

It has been suggested that the radiation of the Crab pulsar is caused by the inverse 
Compton effect. Low energy photons can be scattered by high energy electrons 



An' = T^An 



(3.13) 





(3.14) 
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transferring energy from the electron to the photon. No uniform magnetic field 
is needed to initiate this process. However, the cross section for sufficiently high 
energy transfers is too low to account for the observed brightness temperature 
ruling out this radiation mechanism. 

3.12 Self-absorption 

It is conceivable that the radiation emitted by a plasma is partly reabsorbed. In- 
deed, this so-called "synchrotron self-absorption" is well known [38] and can be 
derived for any source using Einstein coefficients. Below the transition frequency 
which corresponds to the mean particle energy the intensity of the observed spec- 
trum scales as z/^/^ regardless of the particular power law obeyed by the source. 
Because in the radio regime the brightness temperature is many orders of mag- 
nitude bigger than the associated particle energy this effect is irrelevant for the 
understanding of the radio spectrum of a pulsar. 



Chapter 4 

Beam-Beam Interaction in Storage 
Rings^ 

4.1 Introduction 

Colliding particle bunches in a storage ring exert an electromagnetic force on each 
other. The beam-beam parameter ^ is the tune shift exerted by one bunch on 
a particle near the center of the opposing bunch. It is a useful measure of the 
strength of the beam-beam interaction. A limiting value of is reached in an 
e~^e~ collider when further increases in beam intensity lead to particle loss or to an 
increase in the vertical emittance of the beam. In e~^e~ colliders, where the action 
of radiation excitation and damping produces a flat beam, the observed vertical 
beam-beam parameter limit is in the approximate range 0.02 < < 0.1 [121 fT3] . 
At present it is not known whether the emittance increase is due to an incoherent, 
single-particle effect or to a coherent, collective instability of the colliding beams. 
The DCI storage rings at LAL, Orsay, France, used a pair of e"*" and e~ beams 
to collide with another pair, in an attempt to cancel the beam-beam force |14j . 
It was found, however, that the beam-beam limit in DCI was not significantly 
improved by the charge cancellation. Derbenev [15j explained this result in terms 
of a collective instability of the four-beam system and in [ID] the performance of 
DCI was analyzed numerically. This suggests that the beam-beam limit for two- 
beam e~^e~ colliders may also be due to a collective instability. Simulations in 

^This chapter appeared as a journal article [39]. Reprinted in modified form 
with kind permission from the American Physical Society, (c) 2003 by the Amer- 
ican Physical Society 
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|41[ l42l |43] show collective oscillations of the beam at the beam-beam limit. 

In references [151 HBl [E] the stability of the colliding beams was examined 
by solving the Vlasov equation for an equilibrium distribution with small pertur- 
bations. Chao and Ruth [17J considered a beam-beam model in which motion 
was confined to the vertical plane, and in which the beam has a "water-bag" 
equilibrium distribution (uniform within an ellipse in phase space). Synchrotron 
radiation damping and excitation were not considered. When the Vlasov equation 
was solved for a linearized beam-beam force, coherent beam modes were found to 
be unstable near each resonance. In [T6] the stability of a Gaussian equilibrium 
distribution was analyzed with the Vlasov equation for round beams where the 
beam-beam force can be expanded in Bessel functions. A fiat beam model with a 
Gaussian distribution and synchrotron radiation was studied in [1^ |15] under the 
assumption that the distribution always remains Gaussian. A similar approach was 
chosen in |16] for a purely linear beam-beam force. The findings of these models, 
e.g. fiipfiop solutions and period-n solutions are verified numerically in 07| where 
the behavior of fiat and round beams is considered as well. 

In this paper we extend the model of Chao and Ruth to a Gaussian equilibrium 
distribution. In Section 14.21 we set up the equations of motion for the phase space 
distribution and its perturbations, and linearize the beam-beam force. In Section 
I4.3l we solve the equations of motion for radial and angular modes up to first order 
in the displacement from the design trajectory and discuss the implications of our 
results. 
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4.2 Beam Evolution 

We model the flat beam as a current sheet which is uniform in the horizontal 
direction, x, and consider only motion in the vertical direction, y. Consider one- 
dimensional phase space distributions ipiiy, y', s) and 4'2{y, y', s) of the two beams 
which are normalized to unity. Then the deflection from the second (first) beam 
on a particle in the first (second) beam is 

A?/;,2 = -/v2,i(2/,s), (4.1) 

where we define 

I^{y,s) = dysgn{y-y) dy'ilj{y,y' , s) (4.2) 

T J — oo J — oo 

and is the number of particles per unit width in x and is the classical radius 
of the electron. Both beams are assumed to have the same number of particles 
per unit width. The equations describing the motion of iIji^2 are given by the two 
Vlasov equations 

-d^ + y^- ^(^)^^ - (y^ s) = (4.3) 

where the periodic delta function and the unperturbed focusing function are de- 
noted by Sp{s) and K{s), respectively. We want to determine whether the beam is 
stable. That is, we want to know if small perturbations of the phase space density 
grow. Thus, we choose a perturbative ansatz 

ipi,2 = ^0 + AV'1,2 (4.4) 

where ipo is the equihbrium distribution, i.e. a solution of Eq. (14. 3 p with ipi {y, y', s) = 
"02(2/5 y', s) = ipo{y, y', s) = ipoiy, y\ + C), where the circumference of the ring is 
denoted by C. Substituting Eq. (14.41) into Eq. (14. 3p . subtracting Eq. (14. 3 p written 
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for the equilibrium distribution, and neglecting the term which contains a product 
of two perturbations we find 

+ y -Q-y Q^ny. s) - 5.(3)^/a,., = 0, (4.5) 

where 

F{y,s) = K{s)y + 5,{s)^,{y). (4.6) 
If we approximate the beam-beam force as linear in y 

Fiy, s) ^ F{s)y = K{s)y + 5^(^)4 (4.7) 

with 

■ y (4.8) 



y=0 

we can replace K{s) by the perturbed focusing function F{s) to compute the 
perturbed Twiss parameters. In the next step we transform Eq. (14. 5 p to action- 
angle coordinates 

a/2/5 J cos = -./^(sin0 + acos0) . (4.9) 



y 

The betatron function is perturbed by the linearized beam-beam kick from ipQ. We 
form the linear combinations for the a- and the vr-mode 

/± = At/-! ± (4.10) 

Then Eq. (14. 5 p can be decoupled and rewritten in action-angle coordinates as 

The quantity ^ = — ^sin^^-j/'o + ^YJ^'i^'^o^ simplifies since the lineariza- 
tion of the beam-beam force in Eq. (14. 7p leads to ipQ = ipo{J) and we are left 
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with 

df± , ldf± r—— . dipo 1 

^ + ± v^sm05,(.)^/,^ = 0. (4.12) 

In the following discussion we omit the label ±. 

4.3 Solving the Equations of Motion 

When the interaction term in Eq. (14. 3 p is not considered, any different iable distri- 
bution which depends solely on J is an equilibrium distribution. In general, ipo will 
be a function of both J and (p. Fortunately, an arbitrary differentiable function of 
J is an equilibrium distribution, at least to linear order in y after introducing the 
perturbed betatron function. We choose a Gaussian equihbrium distribution 

MJ) = ^^e-i (4.13) 

since in the presence of damping and quantum excitation the beam distribution 
naturally tends to a Gaussian distribution. The deflection of a particle due to the 
presence of a Gaussian beam can be obtained from Eq. (14.21) . 

/*M = ^^erf(^V (4.14) 



7 Vv^ 

We expand the linearized version of Eq. fl4.12p using the ansatz 

OO OO y j\ 

f{J,<P,s) = J2 11 9n'i'{s)e-iL^, (7) (^-^^^ 

n'=OZ'=-oo ^ ^ 

Since the perturbation must be periodic in we can express the - dependence in 
terms of a Fourier series. The orthogonality relation for the Laguerre polynomials 
comes with the convenient weight factor e~i which simplifies working with expres- 
sions that contain the Gaussian equilibrium distribution. Furthermore, using the 
weight factor in the set of basis functions, guarantees that the perturbation falls 
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off as J — > oo. We will refer to the modes represented by the first and second 
index in gni as "radial" modes and "angular" modes, respectively, i.e. these words 
refer to the two-dimensional phase space described by action-angle variables. With 
Eq. (I4.15P the linearization in Eq. (14. 8 p leads to 



If oc 



dy' 



dyf{y,y',s) - / dyf{y,y',s) 

Jy 



(4.16) 



/CO r pO poo 

dy' / dyf{y,y\s)- dyf{y,y\s) + 2yf{Q,i,s) (4.17) 

Using e~^Ln{x)dx = 6no the first part of Ij is given by 

/ dJ d<P[fiJ,4> + 7^,s)-fiJ,<P,s)] = -4ey2go^2i'+i)^jr^. (4.18) 

Jo J-n/2 .trt. + ^ 



r"00 /"""/S 

The second part is given by 



/oo /"OO 1 

#'/(0,y',s) = 2 / dJ-=[f(J,i,/2,s) + f{J,-7r/2,a)] = 



^ n'=0/' 



19) 



Here we have made use of 



OO 

— Lfii^x^dx 





(4.20) 



(2"n!)2 

Inserting into Eq. ( 14. 12^ . projecting this equation onto our chosen set of basis 
functions by means of the orthogonality relation of the Laguerre polynomials 



POO 

/ e~''Ln{x)Lm{x)dx = 6n 
Jo 



(4.21) 



and using 



^ ' 2(2n- l)(2"n!)2 2(2n - 1) 



(4.22) 
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and 

/■oo 

/ xe~^Lnix)dx = SnO - 5nl (4.23) 
Jo 

we obtain 



^ + 'jgni = tMs)^ E f; M^,^n'9n'i' , e = ^a/^^ , (4.24) 
^ ^ n'=oi'=-oo 7 V vre 



where 



1 I'-i 1 ii_ 

-PniSl^i - 5l-i)5n',oi-l)~T:ai> + {5n,0 - ^n,l)('5«,2 " 5/ -2) ^'n' ( "1 ) ~ 

2n — I I' 

= {2'Kl)-^Mnl,n'V (4.25) 

The coefficients ai are 1 for odd I and for even / and vice versa for the 
coefficients hi. Each column and each row of the matrix M refers to one particular 
combination of an n and an / value. 

4.4 Dynamic Tune 



We calculate the tune v in terms of the unperturbed tune by means of Eq. fl4.26p . 

u~uo = ^ i) (3{s){F{s)-K{s))ds . (4.26) 



In order to obtain F{s) — K{s) the deflection in Eq. fl4.14p is linearized. This gives 



._.„.^y!^.5, ,4,T) 

7 V vre 



where /3* denotes the beta function at the interaction point. 



4.5 Coherent Beam-Beam Instability 

We solve the ODE fl4.24p and rewrite the solution in matrix form such that the 
beam transport after one turn is described by a matrix T which acts on a column 



44 

vector G that contains all gnu i.e. G{C) = TG{0). We parametrize the beam- 
current by the linear tune shift parameter ^. One obtains the following relation 
for the Qni^s immediately before and immediately after the interaction point by 
integrating through the interaction point: 

G(0+) - G{0') = ±^MG{0-) (4.28) 

There is no coupling among different Fourier components between collisions. In 
this case Eq. f l4.24l) simplifies to 

which is solved by 

The one-turn transfer matrix becomes 

T± = R{l±iM) (4.31) 

where i? is a diagonal matrix which has the elements e"^'^*''^ on its diagonal. The 
matrix M has the following properties which follow immediately from Eq. (14.251) . 

Mni,n'i' = for / + /' = odd 

Mnl,n'-l' = Mnl^n'l' 
Mn-l^n'V = —Mnl,n'V 

M:i,n'i' = -Mniyi' (4.32) 

In order to decide whether the system is stable or not we have to find out what 
happens to an arbitrary initial perturbation after a large number of turns, i.e. one 
needs to consider the limit where N — > oo. Every matrix norm of the latter 
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quantity tends to infinity if tlie absolute value of one eigenvalue of T is bigger than 
1. To analyze the stability for a given tune v and a beam-beam parameter ^, we 
therefore compute the eigenvalue \max that has the largest modulus. In case of 
instability we compute the corresponding eigenvector G and find its component g^i 
which has the largest modulus. This indicates that the instability mainly drives 
the radial mode n and angular mode /, causing / to be dominated by L„(^)e*''^. 
Since the perturbation / must be real taking its complex conjugate must leave 
/ invariant which gives the constraint Qni = Qn-i- Indeed Eq. ( I4.24p is invariant 
under complex conjugation and replacing / — > —I. It follows that the coefficients 
of T have the property Tn-i^n'-v = Tniyv, which also follows from Eq. (14. 32^ . This 
requires that eigenvalues of T are either real or come in a pair with their complex 
conjugate: Let be a matrix performing the transformation / — > —I then we 
have STSSG = XSG and finally T{SG*) = X*{SG*). Therefore, the /-mode and 
the — / mode are always excited simultaneously with equal strength. 

4.6 Results and Discussion 

In Fig. 14.11 and 14.21 we varied the tune u between and 1 and the beam-beam 
parameter ^ between and 0.12. A point has been plotted if the absolute value 
of all eigenvalues of T is smaller than or equal to 1 for both the a- and the tt- 
mode. We truncated T to the indicated modes. In Fig. 14.11 only the 5 modes 
/ = —2 ... 2 for n = were considered. In Fig. 14.21 we included the same angular 
modes for n = ... 2. The first and second order resonances can be recognized 
clearly. Resonances of orders higher than 2 cannot be expected in our linearized 
model. It is interesting to note that the inclusion of radial modes stabilizes the 
motion of the beam so that a larger ^ can be tolerated. 
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Figure 4.1: Stability diagram for n = 0, I = —2 ... 2 



In Fig. 14.31 and 14.41 we again varied u and and plotted the largest eigenvalue 
|Amax| vs. z/ and determined which mode becomes unstable by selecting the biggest 
component of the eigenvector which is associated with the largest eigenvalue. The 
plot shows that in the absence of dynamics in the radial direction / = ±1 and 
/ = ±2 modes become unstable in the vicinity of u = 0.5, but in Fig. 14.41 only 
/ = ±1 modes are excited around u = 0.5. Furthermore, the unstable / = ±2 
modes which accumulate in the vicinity of z/ = 0.25 and u = 0.75 are attenuated 
if the n = 1 mode is included. Therefore, the radial motion leads to a damping of 
the / = ±2 modes. 

In Fig. 14.51 we computed the phase of the largest eigenvalue of I = ±2 instabili- 
ties, corresponding to quadrupole oscillations (vr-mode only), versus the perturbed 
tune for various Au. The slope of the two lower lines is 2 which indicates that the 
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Figure 4.2: Stability diagram for n = . . .2, I = —2 ... 2 



collective oscillation frequency of the quadrupole mode is twice the single particle 
oscillation frequency for small ^. The spread of the points for fixed u shows how 
strongly the beam-beam parameter ^ influences the frequency of quadrupole oscil- 
lations. In Fig. 14.61 this spread is significantly lower which again shows that radial 
modes have a stabilizing effect. 

The dependence of this spread on u can be understood analytically. For sim- 
plicity we consider only the n = modes. Close to a resonance where lu is integer, 
goi and g^^i perturb the beams the most. Thus, we content ourselves with the 
following 2x2 matrix [17] 
/ 



T 



-2iTilA 



\ 















i 1 


1^ 






1 ± 








I 











(4.33) 



which satisfies all properties listed in Eq. (14.321) for ia = S^Mqi^i. The imaginary 



48 



2 - 



1 . 6 



1.4 - 



1.2 



































■ • 






• • 


• 


• 


• 


• 






■ 


• 


• 


• 






• 


• 


• i 






• • • 


*• • •• 


• 


• 


•• • •* 






• 




*• • •• 


• 






• 


• • • •* 








- • • 


• • • * 
• • 






• • • • 
• • 










. • • 










► • • • 
• • • • 












































• • 


• • 


• 


• 




• 




• 

••* 


• 




• 






• 




•a 
••• •« 
• 


































.. • • 










• • 
















• • • 


•• 






•• 


.* * • 












* m 


• * 


' • • 




• 


■ * 




• • 










' • • 


















• 

••• * 






• 




••I 
















• • 


. ... 


> 


•• < 


... . 


• 








. ... . 








*• * 


>. ... . 










• • • 






• • • 




1 • 






• • ■ 










• • • 










• 






• 




•• 






• 


•.; 






;.• 


• 








• • 


• •* 
• ••• 






••• • 


• •••*< 

• • 


!•* 
• 




• 


• • 
• ••• 










• 

••• • 


• ••••« 
• • 




• 


•••••• • 


• • 

• < 






• • 




m • 


::•.: 




• • 

• 


>»*'• 


• 


• 


'>;> 


• • 

> • 






• 
• 


• ••• • 


• • 


••: 


:•• 


• • 




• 






• • 


m* 


• 
• 


• 

• 


*m 


• • 










* ••• 


■• *• 




••• * 


■ * * 




• 




* ••• 




>• 








• * * 




• 


• 


• • < 






* • • 




• • 






• • 


• 


• 


• 


• 








•• 


• • 


• 


•.: 




• 


• • 


> •• 




• • 


• 




>• 


• < 




• 


• • • 




• 




• • 


• 


• 


■ • 




• 


• 




• • 




• 


• 




• • 








• ••• • 




•••• 


•••• 










*• • 


















• 






1 ■ 


• 4 






• 


• 








• 


• 










• 














• • 










>• 
• 


• 






• • 




•• 


• • * 


• • 






• • 






***'• 


• 


• • 




*• 


•< 




• • 






• 


• • 




• 


• 




■ • 


• 




* • 






• 


• 






• • 




•• 


• • 




1 ••• 


••• * 




• • * 






• • 




* •* 


*• 


•* 


*• * 










• • • 




■ 


• 
















• 


• 












• • 
















• 
































* 


• 








*• 
























■ 




• 


• 














• 




• 


• 


• •••• 


• • 


• • 


••••• 


• 


*.• 




• 


••••• 


• 


• 


• 


• 


••••• 






• 




• 










• 


• 














• 











































0.2 



0.4 



. 6 



. 



Figure 4.3: Absolute value of the largest eigenvalue Xmax vs. tune. Gray points 
indicate unstable / = ±1 modes and black points indicate unstable / = ±2 modes. 
The following modes were included: n = 0, I = —2 ... 2 
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Figure 4.4: Same as Fig. 14.31 but for n = ... 1, / = —2 ... 2 
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Figure 4.5: Phase vs. perturbed tune for = 0, / = ±2 modes (vr-mode only). 

parts of the eigenvalues of the matrix T vanish for eigenvalues whose absolute value 
is bigger than 1. This leads to the plateaus at and 0.5 in Fig. 14.51 and 14.61 at 
tunes V where the / = ±2 mode becomes unstable in the Fig. 14.31 and 14.41 The 
difference between the dipole oscillation frequencies i/^ plotted in Fig. 14.71 gray and 

plotted black of the vr and the a mode divided by the beam-beam parameter 
^ is referred to as the Meller factor |18] or the Yokoya factor |19]. This factor is 
plotted for all points of our computation for which both the vr and the a mode 
indicate stable motion. In Fig. 14.81 one can see that this factor is always above 
1.25 in our Gaussian flat beam model. 

There are only a few points close io v = 0.25 and v = 0.75 since the / = 2 
modes for these tunes are unstable for small ^. 
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Figure 4.6: Phase vs. perturbed tune for n = . . . 1, / = ±2 modes (yr-mode only). 
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Figure 4.7: The dipole oscillation frequencies are plotted gray for the /+ distri- 
bution and black for the /_ distribution with ^ = to 0.2 for n = 0, / = ±1 
modes. 
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Figure 4.8: The Meller factor for stable motion in the region ^ = to 0.2 for n = 0, 
I — ±1 modes. 
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4.7 Possible Extensions 

4.7.1 Higher Order Resonances 

In order to study resonances of order higher than 2 Eq. (14.21) must not be hnearized, 
but rather the double integral has to be expanded about y = to orders higher 
than 1. The expansion to 2nd order contains y"^ j'^^dy' -^f {y^y' , s)\y=Q. Inserting 
the expansion in Eq. (14.151) for / and writing ^ in terms of J and allows the 
evaluation of the integral. The resulting term ^ry^ = —^^{2pj)^e~i sin cos^ 
in Eq. (14.51) needs to be expanded in Laguerre polynomials and gives rise to higher 
orders in radial modes. The n-th order term can be written in terms of powers of 
cos n(f), sin n0 and lower frequency parts. Since the beam-beam force acts only 
at a single point, its contribution is not averaged out in the limit of a large number 
of turns if the tune matches the frequency of one of the sine or cosine functions. 
This is the case if the tune is a rational number, so higher order resonances would 
appear in Fig. 14. 2[ Without truncating the series the model would result in an 
infinite number of resonances since one can always find a rational number between 
two irrational numbers. However, this procedure is complicated by the fact that 
Eq. (I4.13P is not an equilibrium distribution anymore when nonlinear terms are 
included. 

When the length of the bunch and its longitudinal motion is included, synchro- 
betatron resonances can occur [50j when the bunch length is in the order of the 
betatron function. Including these resonances would require and extension of our 
treatment from two to four dimensional phase space. This would be a worthwhile 
but tedious continuation of our work. 
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4.7.2 Damping by Synchrotron Radiation 

One can extend the presented model to account for damping by synchrotron ra- 
diation. In order to obtain the equihbrium distribution in Eq. fl4.13p quantum 
excitation must be included as well. This turns Eq. fl4.3p into the Fokker-Planck 
equation (14.341) . In preliminary computations we found that the graphs we pre- 
sented above remain unchanged for realistic values of the damping and excitation 
coefficients. To simplify the Fokker-Planck equation, we averaged over the phases 
in the damping and excitation terms but not in the beam-beam interaction term. 
This can be justified since the betatron phases in the terms for damping and quan- 
tum excitation change during one turn while the phase in the interaction term 
changes only once per turn. In Eq. (14.341) A is the energy loss per turn due to 
synchrotron radiation divided by the energy of the particle, rj is the dispersion and 
D is the quantum excitation coefficient. 

+ y -t; {T^y + K[s)y + 5p[s)I^^^^ [y, s) 



ds ' dy VC ' ' 1 "^^"^^'^^•^^-'"7 dy' 

-^^.2 + I^(.| + V^)%.2 (4.34) 



4.7.3 Different Tunes 

If the two beams have different tunes, Eq. (14. lOp cannot be used anymore to 
decouple the system. It is easier to work with the uncoupled system and solve 
for the Qni of the two beams separately. Introducing the column vector G which 
contains the gni for both beams, one can proceed as before and describe the beam 
transport for each turn by a matrix multiplication with a matrix T. Introducing 



R 



' R{ui) ^ 
^ R{U2) J 



(4.35) 
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where R{i') is a diagonal matrix which has the components e ^^^''^ we can write 
the matrix T as 

/ 



M 
M 



(4.36) 



Chapter 5 

Coherent Synchrotron Radiation^ 

5.1 Introduction 

The high brightness temperatures of the radio emission of pulsars (Tb ^ lO^^K) 
imphes a coherent emission mechanism [52l [531 EH EH [55] and some part of the 
radio emission of extragalactic jets may be coherent ^56j. Recently, coherent syn- 
chrotron radiation (CSR) has been observed in bunch compressors [HI IIH] 
which are a crucial part of future particle accelerators. When a relativistic beam 
of electrons interacts with its own synchrotron radiation the beam may become 
modulated. If the wavelength of the modulation is less than the wavelength of the 
emitted radiation, a linear instability may occur which leads to exponential growth 
of the modulation amplitude. The coherent synchrotron instability of relativistic 
electron rings and beams has been investigated theoretically by [2n [22l [23l [571 [58] . 
Goldreich and Keeley analyzed the stability of a ring of monoenergetic relativistic 
electrons which were assumed to move on a circle of fixed radius. Electrons of the 
ring gain or lose energy owing to the tangential electromagnetic force and at the 
same time generate the electromagnetic field. [20] analyzed the stability of a rela- 
tivistic electron ring enclosed by a conducting beam pipe in an external betatron 
magnetic field. A distribution function with a spread in the canonical momentum 
was chosen for their analysis. For simplicity the effect of the betatron oscillations 
was not included in their treatment. They find a resistive wall instability and a 
negative mass instability. Furthermore, they find an instability which can perturb 

^This chapter appeared as a journal article [51]. Reprinted in modified form 
with kind permission from the American Physical Society, (c) 2005 by the Amer- 
ican Physical Society 
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the surface of the beam. [57] analyzed the stabihty of a ring of relativistic electrons 
in free space including a small energy spread which gives a range of radii such that 
particles on the inner orbits can pass particles on outer orbits. [58j has developed 
a similar model which includes the effects of the conducting beam pipe. Numer- 
ical simulations by jSH] show the burst-like nature of the coherent synchrotron 
radiation. 

The present work analyzes the linear stability of a cylindrical, collisionless, 
relativistic electron (or positron) layer or E-layer [60]. Particle densities in pulsar 
magnetospheres are very low, of order the Goldreich- Julian charge density ncj = 
0-B/27rce ~ 10^^ cm~''^{B/m^'^ G){R/rY[P{sec)]-^ at radius r > R, where R 
is the stellar radius, B = 10^^i?i2 G is the surface field strength, and P is the 
rotational period; thus, the magnetospheric plasma is collisionless to an excellent 
approximation p8]. The particles in the layer have a finite 'temperature' and 
thus a range of radii so that the limitation of the Goldreich and Keeley model is 
overcome. Although we allow a spread in energies, we assume that it is small, 
so the charge layer is also thin; efficient radiation losses are probably sufficient to 
maintain rather low energy spreads in a pulsar magnetosphere, although the precise 
size of the spread is still not entirely certain. Viewed from a moving frame the 
E-layer is a rotating beam. The system is sufficiently simple that it is relevant to 
electron flows in pulsar magnetospheres (cf. [SI])- The analysis involves solving the 
relativistic Vlasov equation using the full set of Maxwell's equations and computing 
the saturation amplitude due to trapping. The latter allows us to calculate the 
energy loss due to coherent radiation. 

In §5.21 we describe the considered Vlasov equilibria. The flrst type of equilib- 
rium (a) is formed by electrons (or positrons) moving perpendicular to a uniform 
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magnetic field in the z— direction so as to form a tliin cylindrical layer referred to 
as an E-layer. The second type of equilibrium (b) is formed by electrons moving 
almost parallel to an external toroidal magnetic field and also forming a cylin- 
drical layer. §5.31 describes the method of solving the linearized Vlasov equation 
which involves integrating the perturbation force along the unperturbed orbits of 
the equilibrium. In §5.4[ we derive the dispersion relation for linear perturbations 
for the case of a radially thin E-layer and zero wavenumber in the axial direction, 
kz = 0- We find that there is in general a short wavelength instability. In §5.51 
we analyze the nonlinear saturation of the wave growth due to trapping of the 
electrons in the potential wells of the wave. This saturation allows the calculation 
of the actual spectrum of coherent synchrotron radiation. In §5.61 we derive the 
dispersion relation for linear perturbations of a thin E-layer including a finite axial 
wavenumber. The linear growth is found to occur only for small values of the axial 
wavenumber. The nonlinear saturation due to trapping is similar to that for the 
case where kz = 0. In §5.71 we consider the effect of the thickness of the layer 
more thoroughly and include the betatron oscillations. §5.81 discusses the appar- 
ent brightness temperatures for the saturated coherent synchrotron emission. §5.91 
discusses some implications on particle accelerator physics. §5.101 gives conclusions 
of this work. 

5.2 Equilibrium Configuration 
5.2.1 Configuration a 

We first discuss the Vlasov equilibrium for an axisymmetric, long, thin cylindrical 
layer of relativistic electrons where the electron motion is almost perpendicular to 
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the magnetic field. This is shown in Fig. 15.11 The case where the electron motion 
is almost parallel to the magnetic field is discussed below. The equilibrium has 
d/dt = 0, d/d(j) = 0, and d/dz = 0. The configuration is close to the non-neutral 
Astron E-layer of [60]. The equilibrium distribution function /° can be taken to be 
an arbitrary non-negative function of the constants of motion, the Hamiltonian, 

H^{ml+pl+pl + plY^'-e^'ir) , (5.1) 

and the canonical angular momentum, 

= rp^ - erA^{r) , (5.2) 

where = + is the total (external plus self) vector potential, ^'^ is the self 
electrostatic potential, mg is the electron rest mass, — e is its charge, and the units 
are such that c = 1. Here, the external magnetic field is assumed to be uniform, 
= Biz, with Al = rBl/2, and El > 0. Thus we have f = f{H,P^). We 
consider the distribution function 

f = KS{P^ - Po) exp [ - H/T] , (5.3) 

where K, Pq, and T are constants (see for example [62|). The temperature T in 
energy units is assumed sufficiently small that the fractional radial thickness of 
the layer is small compared with unity. Note that a Lorentz transformation in the 
z— direction gives a rotating electron beam. 
The equations for the self-fields are 

i (;^) = I "* • 




Figure 5.1: Geometry of relativistic E-layer for the case of a uniform external 
axial magnetic field. 
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where v^^ = {P^/r + eA^)/H. 

Owing to the small radial thickness of the layer, we can expand radially near 

To 



r 



+ (ro) 



1 

+ 5rDi+-6r^D2, 



(5.6) 



where Di, D2 are the derivatives evaluated at tq, and br = r— tq with (5r/ro)^ <^ 1. 
We choose Tq so as to eliminate the term linear in Sr. Thus, 



where ujpr is the radial betatron frequency, and 



(5.7) 



me < 1 + 



^ + eA^(ro) 
'"0 



1/2 



(5.8) 



Ho 

70 = — 

me 



1 



^^0 



+ eA<^(ro) 



We assume 7q ^ 1 and i^^o > so that i^^o = 1 — l/(27o) to a good approximation. 
The "median radius" Tq is determined by the condition 



2Ho 



dr 



, 



ro 



or 



1 fP. 



Ho \ro 
To a good approximation, 

ro = 



^ + e 



dr 



(5.9) 



rriejov^o 



(1 - 2C)eS 



1 + 2C) or 



2^ 

e5j 



[1 + 3C + 0(C')] 



(5.10) 
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Here, 

C^-^, with C^«l, (5.11) 

z 

is the field-reversal parameter of Christofilos. For a radially thin E-layer of axial 
length L consisting of a total number of electrons A^, the surface density of electrons 
is cr = Nji^Tiri^V) and the surface current density is —ev^a. Because -Bf(ro) 
is one-half the full change of the self-magnetic field across the layer, we have 
C = reN/{'^L), where re = e^/{mc^) is the classical electron radius. Notice that 
iV, and 7L are invariants under a Lorentz transform in the z— direction. 
The radial betatron frequency uj^r is given by 

Hi 1,.^^ ■ (5.12) 



Using Eq. (15.91) gives 



2 .2. g ^2^s 



1 - 2C C 



(5.13) 



Tg roAr7^ 

The term oc 1/Ar is the sum of the defocusing self-electric force and the smaller 
focusing self-magnetic force. For the layer to be radially confined we need to have 



C < A/vr/2 7^(Ar/ro). For ( <^ 7^(Ar/ro) and ^ 1, we have cu/j^ = 1/ro to a 
good approximation. 

The number density follows from Eq. (15. 3p . 

/ \ 1/2 

n ~ no exp — - where Ar 



2Ar2 y \ 



2C- v^C(ro/A077^ 



or ^ (5.14) 
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where 



rp \ 1/2 



vth = (5.15) 

and 

no = 27r KHoTr^ exp I I . (5.16) 

As mentioned we assume the layer to be radially thin with (Ar/ro)^ <^ 1. Conse- 
quently, Eqs. (15. 4p and (15.51) become 



AnenQ exp 



^^2 47reno t;</,o exp (^-^J . (5.17) 

Thus we obtain 

C - — si ^^-^^^ 

The equilibrium is thus seen to be determined by three parameters, 

C\ vl, and l/7o', (5.19) 

which are all small compared with unity. 

5.2.2 Equilibrium Orbits 

From the Hamiltonian of Eq. (15.71) we have 
(P5r 



ujfi^6r, 6r{t') = sm[uji3r{t' - t) + ip] , (5.20) 



where r — ro = Svi sin ip. For future use we express the orbit so that r(t' = t) = r 
where (r, t) is the point of observation. Also, we have 
d(h Ps + erAJr 



dt ruejr'^ ~^ dr 



5r+.. , (5.21) 

ro 
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so that 



0(t') = + (t' _ + 



UJpr Or 

— 5ri coii[uj prit' — t) + (p\+ 5ri cos(v?) 



(5.22) 



where 90/9r|ro = —(Pq/tq. For ^ <^ 7^(Ar/ro) and ^ 1, we have d<^ldr\r^lujfir = 
— I/tq to a good approximation. Because the E-layer is uniform in the z— direction, 



z{t') = z + {t' - t)v, . 



(5.23) 



The orbits are necessary for the stabihty analysis. 



5.2.3 Configuration b 

Here, we describe a Vlasov equihbrium for an axisymmetric, long, thin cylindrical 
layer of relativistic electrons where the electron motion is almost parallel to the 
magnetic field. The equilibrium distribution function /° is again taken to be 
given by Eq. (15. 3p in terms of the Hamiltonian, if, and the canonical angular 
momentum, = rp^ — erA^{r), where = A^. We make the same assumptions 
as above, 7^ ^ 1, T/^nie'y) <^ 1, and Ar^/r^ <^ 1. In this case there is no 
external field. Instead, we include an external toroidal magnetic field with 
corresponding vector potential Al and an external electric field with potential 
The fields B*^ and correspond to the magnetic and electric fields of a distant, 
charged, current-carrying flow along the axis. Thus, \Ef.\ < \B^\. The considered 
external field is of course just one of a variety of fields which give electron motion 
almost parallel with the magnetic field. Note also that the distribution function 
is restricted in the respect that it does not include a dependence on the canonical 
momentum in the ^— direction = rrie'jVz — eA^. 
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The distribution function (]5.3p gives = so that there is no toroidal self 
magnetic field. Thus the self-potentials in this case are also given by Eqs. fl5.4l) 
and (15.51) . Eqs. (15. 6p - (15.91) are also applicable with the replacement of by the 
total potential $. In place of Eq. (I5.10p we find 

- (1 - 20emro) ~ i^^' + ' ^'-'^^ 

where ( = BKr^)/ E^.{rQ). We again have ( = r^N/ (7-L), where = e^/ (mc^) is the 
classical electron radius and L is the axial length of the layer. Because /dr"^ = 
— (l/r)(i$^/(ir, the radial betatron frequency is again given by Eq. ( 15.13^ (with $ 
now the total potential) so that the equilibrium orbits given in §5.2.21 also apply 
in this case. The electron motion is almost parallel to the magnetic field in that 
{Bl/Blf = C^{E^/BlY < < 1. Notice that Eq. (15:241) for tq is formal in the 
respect that E!^ oc 1/r. Therefore, tq is in fact arbitrary in this case. Because the 
wavelengths of the unstable modes are found to be small compared with tq, it may 
be interpreted as local radius of curvature of the magnetic field. 

5.3 Linear Perturbation 

We now consider a general perturbation of the Vlasov equation with /(r,p,t) = 
/"(r, p) + 6f{r, p, t). To first order in the perturbation amplitude 6f obeys 

where 5E and 5B are the perturbations in the electric and magnetic fields. All 
scalar perturbation quantities are considered to have the dependencies 

F(r) exp(im0 + ikzZ — iut) , (5.26) 
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Figure 5.2: Geometry of relativistic E-layer for the case of an external toroidal 
magnetic field with an external radial electric field. 
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where the angular frequency u is taken to have at least a small positive imaginary 
part which corresponds to a growing perturbation. This allows for a correct initial 
value treatment of the problem f63]. For a perturbation taken to vanish as t — > — cxd, 

dp 



5f{r,p,t) = e J dt'l^6E[r{t'),t']+v{t')x6B[r{t'),t' 



dp 



(5.27) 



where the integration follows the orbit [r(t'), p{t')] which passes through the phase- 
space point [r, p] at time t. For the considered axisymmetric equilibria, 



d 



^ - P. 
'dp ~H 

p d p dPs d 



dp HdH H dH dP^ 



df p dP^ / df 



since f = f{H, 



dp HdH \dP^ 
and H = H{P^, . . .) 



+ 



H 



dH df 
dK dH 



(5.28) 
(5.29) 

(5.30) 



a d 

V 



dp dH 
df p df 



d 



dp H dH 



dP, 
df 



H 



dPs 



(5.31) 
(5.32) 



H 



where the partial derivatives are to be evaluated at constant P<^ and H, respectively. 
Thus, the right-hand side of Eq. (I5.25p becomes 



d5^ ■ \ df 

+ iujU 6^1 - 5$) + iuj\±_ ■ 6 A + 

/ dH 

df 



dt 

d5^! 



dt 



+ im{(j) 5^ — (5$) + imv_L ■ 6 A 



dPs 



(5.33) 



where 5E = — V5$ — d6A/dt and 5B = Vx5A, 5\E' = r6A^ is the perturbation in 
the flux function, v_|_ = {vr, v^), and d/dt = d/dt + \ ■ V. We assume the Lorentz 
gauge V-5A + d5^/dt = 0. 
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Evaluating Eq. (I5.27P gives 



+ e 



dH 



oo 

t 



-5$ + lu / dt' 6' 5-^' - 5$' + v'l ■ 5 A' 



+ / dt' UW-5<l>' + v'i -5 A' 



(5.34) 
(5.35) 



where the prime indicates evaluation at \r{t'),t']. The integration is along the 
unperturbed particle orbit so that dp /OH and df^/dP^ are constants and can be 
taken outside the integrals. Note also that d/dt acting on a function of (r, t) is the 
same as D/Dt. 



5.4 First Approximation 

As a starting approximation we neglect (i) the radial oscillations in the orbits 
[(Ar/ro)^ ^ 1], (ii) the self-field corrections to orbits proportional to C, (hi) the 
terms in 5f proportional to (f^\ ~ (Ar/rg)^ ^ 1), (iv) we take fc^ = and (v) 
we assume the layer is very thin. Owing to approximation (iii), we can neglect the 
terms oc v_l -(JA in Eq. (15.351) in the evaluation of 5p and 5J^. This is because these 
terms give contributions to 5f which are odd functions of Vr and Vz- Therefore, 
their average contribution can be neglected. 
Evaluation of Eq. (15.351) gives 



5f = -e 



dp 



dH 



P4, ^-m 



dPs 



(5.36) 



where (f) = <p{ro)- The approximations lead to a closed system with potentials 
((5$, (5\E') and sources {6p, 5J^). 
We have 

5p= -e j d^p 5f = j dprdpzdP^y Sf , 
5J^ = ~^ J ^^P ^0'^/ = ^~ J dprdpzdP^ v^6f . (5.37) 
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For the considered distribution function, Eq. ( 15. 3p . df^/dH = —f^/T. The 
dp /dP^ term in Eq. (15.361) can be integrated by parts. Furthermore, note that 
dH/dP^ = <p and d(f)/dP^ = — (0)^/if, which corresponds to an effective "negative 
mass" for the particle's azimuthal motion [6ll[6l[65]. From the partial integration 
the small term proportional to dv^/dP^ = v^/Ij-qH^) is neglected. Also note that 
H is not a constant when performing the integration over momenta. Evaluating 
this term by an integration by parts with a general function g{P^) in the integrand 
gives 



dP^ 



dP^ 



H 



- K j dP^5{P^ - Po)^ [9{P<^)e-'"^] = 

- K j dP,S{P, - Po)^ [g{P,)] e-"'^ 

+ dP,6{P^ - ^o)^?(P^)e-^/^|^ . (5.38) 

That is, the integration produces an additional term which cancels the 1/T-term. 
Thus, 

j dP^^f = y dP^— ^^-p , (5.39) 

where Aa; = u — m(j). Integrating over the remaining momenta gives 

(5p, 5J4,) = {l,v^) n(r)-^ ■ ^ . (5.40) 

For a radially thin E-layer we may take 



n{r) = no exp {-6ry2Ar') ^ noV27rAr 6{6r). (5.41) 

We comment on this approximation below in more detail when we include the 
radial wavenumber kj. of the perturbation. Then Eqs. (]A.4p and (]A.5|) can be 
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written as 



[5^ro), S^iro)] = [l, rov^il + Au)] 27rVo Z j dr 5p{r) , (5.42) 

where Z = iJm{ujro)Hm{i^rQ), uj = uj/{m^) and Au = Auj/{m(p). Integrating 
Eq. (15.401) over the radial extent of the E-layer and cancehng out the field ampli- 
tudes gives the dispersion relation 

1 = 27rVo hoe^v^Ar Z^ ' . (5.43) 

In terms of dimensionless variables this becomes 

1=.CZ (2A^-^)(^. (5.44) 

where Z = iJm{'mujv^)Hm {mujv^) , Hm = Jm + 'iYm, and the field-reversal param- 



eter ( = AnenQV^Ar ^J-k j^jB^ as given by Eq. (15.111) . 
For m ^ 1 approximation (1B.2P can be used to give 



Jj^mujv^ ^ (2/m)^/^Ai(u;) 
Y^{m^v^) ^ -{2/my^^Bi{w), (5.45) 



where 



w = (m/2)2/3(7"2 - 2Au). (5.46) 



Thus we have 



Z = tJmH^^ ^ (2/m)2/3[Ai(w) Bi{w) + iAi\w)]. (5.47) 
Occasionally, Zm{w) is denoted by Z. For ^ 1 using fIB.SP 



Z ^ {2/mf/y{27T\w\^/^). (5.48) 
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and for \w 



2 < 0.5 using (EH) 



— c^w"^) + i{ci — C2W 



2- 



(5.49) 



For IwP <^ 1 



Z ^ (0.347 + 0.200 i)/m^/^. 



(5.50) 



5.4.1 Range of Validity 

We are interested in the regime where the wavelength of the emitted radiation is 
comparable to the "bunch length", i.e. u; ^ m or equivalently Auj <^ 1. However, 
Eq. f l5.44p is only valid if Alj <^ 7"^. Since we neglected 6Jr and 6Jz we obtain 
from the continuity equation SJ^ = ^Sp. Due to this approximation the factor 
on the right hand side can become bigger than the speed of light if Auj > 7"^ 
which leads to unphysical results. In the latter case SJ^f, = v^6p is a better approx- 
imation. Fortunately, Auj <^ 7~^ is the most interesting case and in the remainder 
of this paper we will always work in this limit. Furthermore, for the continuum 
approximation to be valid the mean particle distance has to be much smaller than 
the wavelength. 

5.4.2 Growth Rates 

It will prove useful to define two characteristic values of m: mi = C^l'^^y'^ and 
777.2 = 27^, and therefore mi = C^^l'^m^jl. We can obtain approximate solutions 
to Eq. fl5.44p in two different cases. There may be solutions with small values of 
7^Aa}, so that w ~ {mjm-})'^!'^ . In this case, Eq. (15.441) becomes a simple quadratic 
equation, which can be solved for Au. We can simplify the solution somewhat by 
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changing variables to a = 7^AtD in which case Eq. (15.440 can be written in the 
form 



1 = 1^(^-1) 



where we have neglected a compared to one in the approximate version of this 
equation. We find that 



a ^ ^/~nCZm7^ . (5.51) 
For case I let us assume that m <^ m2, in which case Eq. (15.511) implies 

a - ±1.121(mi/m)^/=^ e'(Wi2) 
= 1.121(mi/m)^/3(-0.2588 + 0.9659i) . (5.52) 

so |cr| <^ 1 for m 3> mi. The growth rate of the unstable mode is 

l.O83C^/^m2/30 
Ui ^ (5.53) 

7 

in this regime. For case II we assume that m 3> m2, in which case Eq. (15.511) 
implies 



a 



^1/2 

^1/2 

^1/2 



note that the growth rates in cases I and II match almost exactly at m = m2, 
where |cr| ~ C}^'^ . 

Note that is the approximate frequency of the peak of the single particle 
synchrotron radiation spectrum. For more accurate results we employ a numerical 
method for solving Eq. (15.441) outlined in |66|. This method also allows us to count 
the number of roots which are enclosed by a contour. The basic idea is that for a 
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nuU-homotopic cycle T which does not cross any poles or roots and a meromorphic 
function / which is not constant [67] 



where A^(0) is the number of roots minus the number of poles enclosed by F (an 
n-th order root or pole counts as n roots or n poles, respectively). So far we have 
no numerical evidence of the existence of more than one solution with a positive 
real part. The numerical results agree very well with our approximations even if 
m < mi and are shown in Fig. 15. 3[ 

5.4.3 Comparison with Goldreich and Keeley 

Goldreich and Keeley pi] find a radiation instability in a thin ring of relativis- 
tic, monoenergetic, zero temperature electrons constrained to move in a circle 
of fixed radius. Under the condition 1 ^ m}/"^ <^ 7 their growth rate is cUi ^ 
1.160 m^^^lreN / ('j^ro)Y^'^ which is close to our growth rate with L replaced by tq. 

5.5 Nonlinear Saturation 

Clearly the rapid exponential growth of the linear perturbation can continue only 
for a finite time. We analyze this by studying the trapping of electrons in the 
moving potential wells of the perturbation. For (Ar/ro)^ ^ 1, the electron orbits 
can be treated as circular. The equation of motion is 




(5.55) 



e[(5E^ + (vx5B)^] , 



(5.56) 



dt 



where is the canonical angular momentum, where 



eSE^Q exp{ujit) cos(m<^ — Urt) 



(5.57) 
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Figure 5.3: The graph shows the frequency dependence of the growth rate for a 
sample case where 7 = 30 and C = 0.02 obtained from our approximations for 
Eq. fl5.44p . For these parameters, mi ^ 10^ and m2 ~ 2.7 x 10^. 
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where 6E^q is the initial value of the potential, ujr = Re(u;), and uji = lm{uj) 
For a relativistic particle in a circular orbit, 

„3 



SP<p = me*rl54>, where me* = ^^'''^^ ~ —mej, (5.5^ 



where me* is the "effective mass," which is negative, for the azimuthal motion of 
the electron ([MIE] or [65], p. 68). Combining Eqs. (15.561) and (15.581) gives 



-uj^{t) sin (f , (5.59) 



where ip = mcj) — Utt + |7r, cut = luto exp(co'it/2), and cuto = [em(5-E'0o/("^e7'"o)]^''^, 
where ut is termed the "trapping frequency." At the "bottom" of the potential 
well of the wave, simp ^ (p. An electron oscillates about the bottom of the well 
with an angular frequency ~ ut- This is of course a nonlinear effect of the finite 
wave amphtude. A WKBJ solution of Eq. (I5.58p gives 

oc (^2^Q''^exp(— i^jt/4) sin |(2u;To/i-^i)[exp(ci;j)f:/2) — 1]} . (5.60) 

The exponential growth of the linear perturbation will cease at the time tsat when 
the particle is turned around in the potential well. This condition corresponds to 
'^riisat) ~ ^i- Thus, the saturation amplitude is 



/ me7 y / uJijm) 
V erom / 



|5E,,/=(^) (Zli^) , (5.61) 



where \5Esat\ = \SE{tsat)\ = \SEo\exp{ujitsat)- 

5.6 First Approximation with kz ^ 

Here, we consider ^ but keep the other approximations. Our ansatz for Sf is 
general enough to handle this case since it retains the biggest contribution to the 
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Lorentz force in the z-direction which is of the order Vf^Br- In place of Eq. (I5.39P 
we obtain 



j dP^Sf = -eJ dP^ 



(5.62) 



H {uo — mcj) — kzVzY 
where we assume without loss of generality A;^ > and m/ro,uj. In place of 

Eq. fICTll we find 



e{uj,kz) = 1+ kzA{uj,kz) / dv_ 



^ Vth 



[..] = 0, (5.63) 



where 



[...] 



mcf) 



{ijj — mcf) — kzVzY 
Here, e acts as an effective dielectric constant for the E-layer, and 



A{u,kz) =7: C Z{uj,kz) (u- 



kd 



LU — mc 



u = 



(5.64) 



k-zTt^ ) ky 
Z ^ rUr,{u- - klfl-\ if£)[ro(.;^ - A:^^/^ 

and k^ = m/rQ is the azimuthal wavenumber. The expression for Z is from §5.41 
An integration by parts gives 

exp(-f^/2fj\) m(pvjkz 



e{u) = 1 + A{uj) / dv. 



2vr vl 



(5.65) 



where the k^ dependence of e and A is henceforth implicit. We can also write this 
equation as 



eiu) = 1 + B(u) 



u 



Vth 



u 

Vth 



(5.66) 



where 
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and 



/27V J X-Z 

for Im(z) > 0, and 

I r , exp(-a;V2) . ^ f 
V27r i-oo a; - z V 2 

for Im(z) < 0. The second expression for F{z) is the analytic continuation of the 

first expression to \m.{z) < which corresponds to wave damping (see, e.g., |68j . 

ch. 5). Note that terms of order Ao) have been omitted. 

For m ^ 1, the factor Z = iJm{Jm + *^m) can be expressed in terms of Airy 

functions in a way similar to that done in §5.51 One finds Jmlroico"^ — k'^Y^'^] ~ 

(2/m)V3Ai(M;), Y^[ro{uj^ - A;2)i/2] ~ -{2/my/^Bi{w), 

(2 \ 2/3 / 2 \ '^^^ 

— Ai(w)Bi(w) , Zi^ { — ] Ai\w) , (5.68) 
m / \m J 

where 

kz /m\2/3 / I \ 

tan-j/y = — , w = — ^ + tan ^ ~ 2m tan ^/^ 
\2J\Y J 

It is clear that e has in general a rather complicated dependence on -u = -Ur + iui 
and tan?/). Note that the expression for w goes over to our earlier w for = 
noting that utsmip Auj. 

A limit where Eq. (15.661) can be solved analytically is for |mP = | Aa}^/ tan^ ip ^ 
f^^j, that is, for sufficiently small tanip. In this limit Eq. (15.661) can be expanded as 
an asymptotic series F{z) = —1 / z — 1 / z^ — 3 / z^ — ... Keeping just the first three 
terms of the expansion gives 



e 



For tanip and Auj ^ 7 ^, this is the same as Eq. (15.441) as it should be. In 
general Eq. (15.691) will have more than one unstable mode. In the remainder of this 
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0.15 I I , I , I 




-0.1 I I 

0.01 0.02 0.03 0.04 0.05 



tan = / 




tan i|j = / 



Figure 5.4: The figure shows the growth / damping rate uji and real part of the 
frequency Auo,,. = uj — mip in units of as a function of ianip = k^/k^ for m = 100 
and m = 1000 for an E-layers with 7 = 30, ( = 0.02 and Vth = 30/7^. In the 
region of damping Ui < 0, the second expression for F{z) in Eq. fl5.66p is used. 
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dei / dur 


dei/dui 







paragraph we will only study the largest unstable solution for which we recover 
the growth rates found in §5.41 in the limit tan^ — > 0. Fig. 15.41 shows some sample 
solutions. For the case shown the u dependence of Z is negligible. 

General solutions of Eq. (15.661) can be obtained using the Newton-Raphson 
method ([SS], ch. 9) where an initial guess of {ur,Ui) gives (er,ej). This guess is 
incremented by an amount 



(5.70) 



and the process is repeated until Sr = and Ei = 0. Fortunately, the convergence 
is very rapid and gives \e\ < 10^^° after a few iterations. 

Fig. 15.41 shows the dependence of the complex wave frequency on the tangent of 
the propagation angle, tanip = kz/k^, for a sample cases. The maximum growth 
rate is for ip = or k^ = 0. With increasing ip the growth rate decreases, and for 
ip larger than a critical angle ipcr there is damping. For the damping the second 
expression for F in Eq. (I5.66P must be used. Roughly, we find that the critical 
angle corresponds to having the wave phase velocity in the direction of the order 
of the thermal spread in this direction, that is, Ur = Aur/kz ~ Vth- This gives 



(. 



TeN 



1/2 



1 1 

< 



(5.71) 



for nil < 1TL < 1712- Note that the dimensionless parameter which determines 
the cut-off at tsmipcr is I'^Vth- Our numerical calculations of ipcr give a slightly 
faster dependence, ianipcr oc l/m°'^'' for this range of m. Fig. 15.51 shows the m- 
dependence of the critical angle. It is reasonable to assume that in a particle 
accelerator the weak focusing in the z-direction sets a low limit on kz- 




Figure 5.5: Critical angle for 7 = 30, C = 0.02 and Vth = 30/7^. 
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5.7 Nonlinear Saturation for /c^ 7^ 

We generalize the resuhs of §5.61 by including the axial as well as the azimuthal 
motion of the electrons in the wave. The axial equation of motion is 
(Pz 

~ — e^E'^o exp(Li;it) cos(m0 + /c^z — cut) . (5-72) 

The approximation involves neglecting the force oc Vr5B^ which is valid for a 
radially thin layer (Ar^/rg ^ 1). Following the development of §5.61 the azimuthal 
equation of motion is 

nT'el'To—— = —e6Ezo cos(m0 + kzZ — out) . (5.73) 

Combining Eqs. fl5.72p and fl5.73p gives 



d Lf e mSE, 



00 



^„ (1 + tan^ ^) sin , (5.74) 

where (p = m(j) + kzZ — ujt+^TT and tanip = k^/k^. Because <^ 1 for wave growth 
(Eq. fl5.7ip ). the saturation wave amplitude 5Esat is again given by Eq. fl5.6ip . 

5.8 Thick Layers Including Radial Betatron Oscillations 
5.8.1 The Limit /c^Ar > 1 

In this section we include the small but finite radial thickness of the E-layer. We 
keep the other approximations mentioned at the beginning of §5.4[ In particular 
we consider k^ = 0. In order to include the layer's radial thickness, we consider 
the wave equations within the E-layer, 

(V^ + cj2)5$ = -47r5p , 

(V^ + cu^)^^ = -AnrdJ^ , (5.75) 
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where 

is the adjoint Laplacian operator. 

Within the E-layer, we assume that the potentials can be written in a WKBJ 
expansion as 

(5$, 5^') = K<^) exp [im0 + ikr{r - tq) - ioot] , (5.77) 

where kr is the radial wavenumber with (/crAr)^ ^ 1 (fC^,^^) are constants. 
This is equivalent to assuming that the charge density is constant between ro — Ar 
and ro + Ar and zero elsewhere. Evaluation of the time integrals in Eq. fl5.35l) for 

T = Tq gives 



oo 

-iz sin " 



* = J2 e-'"*^n(^) (5.78) 



t 

^jJ ^^^-iujf +im4>' +ikrr' 



oo 



/dt'5^ exp < —itot' + im + (t' — t)0o- 
oo 

+ ikrSvi sm{ujpr{'t' — t)) + ikrf-Q 
rft'5$e"'"*'+*™^°+*''-'^''-*™^ e-^"-^'-*'j„(-A;,5r,) (5.79) 

-OO 



St ■ 

— {- COs[uJprit' - t)] + 1) 



r = 5$(ro,t) V -^n(fc^^-)^"^^P(-^M^--^^^) , (5.80) 

J-oo «(m0 + nujpr - uj) 

where n is an integer, k = (/c^ + k'^Y^'^, with A;,^ = m/ro, and tan-?/' = k^/k^. There 
is an analogous expression for the integral of 5^ . We have used Eq. (15.201) for the 
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radial motion with ip = assuming <^ 1 and ( <^ 7^(Ar/ro) so that = I/tq, 
and Eq. (15.221) for the 0-motion with d(l>o/dr\rg/uj[^r = —^/ro- Using Eqs. (15.351) 
and (I5.80p . the momentum space integrals (15.371) can be done to give 



+ 



H [mcf) — ujY 

Jo(Mr,)-l + (m0-.;) V 

J I — m0 + ncuisr — 



n=—OD 
oo / 



m{K^(f) - K^) ^ — + 2^ — S ,(5.81) 



H [ (m0-a;)2 {m(j) + nu^r - 1^) 

and finally if Auj <^ 7"^ 

e^nQmcjP'Ki^ / Tq^cj — m[l — (1 — Fo)/7^] m ^ ' F„ 



H \ (a; - r7j0)2 ^^„n^ ("^0 + '^^^/Jr - tu)^ 



(5.82) 



The prime on the sums indicate that the n = term is omitted. Here, 



with 

X = kAr . (5.84) 
The 1/T terms in Eq. (I5.82p do not cancel exactly. They may be neglected if 

|A(I;|2 < Fovl^ (5.85) 

for the n = term or if 

\Au;\\n/m - Aoj\ < (5.86) 

for the n 7^ terms. 

For weak E-layers we have for x ^ 0, Fq ^ 1 and Fn^Q — > 0. In this limit we 
recover the results of §5.4[ For x ^ 1 1 ^ ^r'^o ^(/-''^O; the Gaussian factor 
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in the integrand of F„ can be neglected so that one obtains 

Fn ~ I e ^ , --— cos — , even n 

V27rx ^ 2 / 



_ •« sin — , odd n . (5.87) 

An ahernative approximation for F„ can be obtained by using the integral 
representation of the Bessel function. The remaining integral can then be computed 
numerically more easily. In this way we find 

Fr, = / d^exp \-zne - {x^/2){kA,/k - smOf] . (5.88) 

2vr 

For X ^ 1, 1 ^ {k^/krY and \n\ < ^Jx we can approximate sin 6' in the exponent 
by a parabola at its maximum. We obtain 

•?ig— m</) 

Fn ^ . (5.89) 

23/4r(|)yx 

In general F„/ (z"'e~*"'^) decreases as x n increase. This acts to prevent the 
unlimited increase of the growth rate as m — > oo, and it ensures that the sums 
over n converge. Fig. 15.61 shows a plot of Fq obtained by numerical evaluation of 
Eq. dEHH]). 

Within the E-layer, Eq. (15.751) gives 

2 ^ AT^e^nQiTLcj)^ 



k: = uj^ ^ + 



H 

ro0cj — m[l — (1 — Fo)/7^] m F, 



X . 

{uj - m0)2 7"^ ^ (m0 + nujfjr - 00^ 

In terms of dimensionless variables this equation becomes 



kl = 2m^Au - ^ + 



2 



CV2 



7 v^VthV^ 



X 



;i + A^)(l - 7-^) - [1 + (Fo - l)/72] _ ly' j ^ 



(Ac<j)2 72-^ (t>^ n/m — Acij)^ 




Figure 5.6: Fq for Vth — 0.01 and k^ro — 10^ 
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where kr = rokr, = r^k^, k = r^k, and x = kvth- 



Notice that Eq. (15 .77^ can also be written as 



5$ = C2 sin [kr{r - Tq)] + C3 cos [kr{r - Tq)] , 



(5.92) 



for To — Ar < r < tq + Ar. For r < tq — Ar, we have 



5$ = CiJ^(cjr) , 



(5.93) 



since the potential must be well behaved as r — > 0. For r > tq + Ar, we must have 



This combination of Bessel functions gives 5$(r — 00) — for the assumed con- 
ditions where Im(ti;) > 0. Note that these potentials are just the solutions of 
Eq. (]5.75p in our approximation for 6p. The eigenvalue problem can now be solved 
by matching the boundary conditions. However, we have not solved the full eigen- 
value problem. Instead we consider unstable solutions with the restriction that 
krAr ^ 1. Under this condition we can interpret Eq. (]5.9ip as a local dispersion 
relation. Unstable modes found from Eq. (15.911) will need a slight correction in 
order to satisfy the boundary conditions. 

We expect that Eq. (I5.9ip has solutions near each betatron resonance at AO = 
±.71 /m. This is a familiar concept in the treatment of resonances in storage rings 
(cf. [T7] or [39]). We extract each solution by summing over a single value of n 
and —n only and obtain from Eq. (]5.9ip for the case n 7^ and Auj ^ 7^^ 




(5.94) 



7^^i/iV^ [kl + m^'-i 2) 




Thus, 



Auj 




(5.95) 
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for sufficiently big S, i.e. we expect the imaginary part of Au; to be negligible for 
the n 7^ modes. Despite a lot of effort we were not able to prove this statement 
under more relaxed conditions. 

We can easily find an analytic solution of Eq. (15.911) for the case where the 
n = term is dominant. If |Au;| ^1/7^ and |Au;| <^ Fq/'j'^, we obtain 

ACo = ± ^ '^'^^^ (5.96) 

The dependence of the growth rate on kr becomes significant when 7^A;^/m^ 
is comparable to unity. For m ~ m2 ~ 7^, we see that this happens when 
{krArY /'-f^v'ff^ ~ 1, which involves the combination 'y^Vth again. 

The growth rate of Eq. (I5.96P is proportional to y^. This implies from §5.5l that 
the emitted power scales as the square of the number of particles in the E-layer 
which corresponds to coherent radiation. Sample results are shown in Fig. 15. 7[ We 
conclude that the main effect of the betatron oscillations is an indirect one. The 
radial motion itself is unimportant for the interaction. However, the influence of 
the radial motion on the time dependence of the azimuthal angle of a particle is 
important since a shift in can take the particle out of coherence with the wave. 
This effect is accounted for by Fq. 



5.8.2 Qualitative Analysis of the Effect of the Betatron 
Motion 

Let us suppose that Vth ^ 1/7^5 ^^at \Ail!\ is not necessarily small compared 
with Vth (We can still assume \Auj\ <^ 1 without requiring the more restrictive 
condition 7^|AtI'| ^ 1.). The key effect of the betatron oscillations is to "wash 
out" the phase coherence of the response within the layer; for a cold layer, all 
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orbiting particles move in "lock step" , which is particularly favorable for a bunching 
instability. Let us suppose that |Au;| has a real part that is substantially larger 
than 1/7^. The response in the layer scales as an Airy function with argument 
+ where |^| < Vth- The phase accumulated across the layer thickness ~ Vthro 
is 7] ~ ^'^th if ^ '^ih and rj ~ mvlj^ {^Cor / VthY^"^ if AcD,. ^ Vth- Large rj 
ought to imply substantial decoherence of the response in the layer. We see that 
this is likely irrespective of the value of Acj^/^t/i provided that m ^ '^th'^^ i-e. 
for ^ [ri^Vth]^'^^'^ ■ At large values of 'j'^Vth, phase smearing should suffice 

to suppress - if not eliminate - the bunching instability at frequencies near the 
synchrotron peak. Moreover, if ■j'^Vth ^ C^^; the instability should be suppressed 
over the entire range m > for which we found unstable modes in §5.41 

Large Auj^/vth would merely accentuate the smearing. At a given value of m, we 
see that Aur ^ i^^'^'^Jth)"^ ^ i-^- 7^Aa)r > {m/ •y^)^^ {•y'^vth)''^ suffices for large phase 
decoherence in the layer. 

5.8.3 The Limit KAr <C 1 

In order to determine the lowest allowed value for kr and the highest possible 
growth rate the full eigenvalue problem has to be solved. We estimate the result 
by evaluating Eq. (lA.4p in the thin approximation again. Looking at Eq. (lA.4p and 
replacing the Bessel functions by their Airy function approximations for the case 
nil and m <ti m2 we see that the thin approximation is justified if krVth ^ 1 
and m?^^Vth <^ I. It starts to fail completely if m^^^Vth ^ 1, i-e. once we start 
integrating over the oscillating and/or the exponentially damped/increasing part 
of the Airy function, which implies we would like to have m^/^|AtI;| <^ y/Fo with 
|AtDp ^ ^o^th from the previous paragraph. However, for real values of kr we 
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Figure 5.7: Growth rates in the hmit krVth ~> 1 for our reference case 7 = 30, 
C = 0.02 and Vth = 1/7^ and various values of kr- The hne proportional to m ' 
is shown for comparison. 
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expect that the thin approximation will still give us an upper bound of the growth 
rate because it is easier to maintain coherence if all the radiation is emitted from 
the same orbit. With Eq. (15.821) we obtain in the limit Auj <^ 7^^ 

The growth rates can be found as before. For m 3> mi we obtain 

uJi ~ ^ -^/Fo 5.98 

7 

and 

^1/2^1/2^ 



(jji 



7I/2 



for m ^ 7712, i-e. there is an additional factor of v'^. The results for our reference 
case are plotted in Fig. 15.81 which were computed numerically. In Fig. 15.91 the 
function Fq is plotted which we compare with the squared ratio of our new growth 
rates to the ones evaluated previously without betatron oscillations. 

We could also study the effect of the non-zero thickness alone without betatron 
oscillations setting Fq = 1 and -F„^o = and solving the full eigenvalue problem. 
Due to the complicated nature of the dispersion relation we have not done this yet. 
Note that the thin approximation will suppress certain modes, e.g. the negative 
mass instability cannot be expected to be present with the fields having been 
evaluated at one radius only, cf. [7]. 

5.9 Spectrum of Coherent Radiation 

Having computed the growth rate and the saturation amplitude, the radiated power 
can now be calculated. Starting from Eq. flA.lOl) we now have 



Pm — ij^Lur'^ 1 57^0 1 



2 



(5.99) 
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10^ 10^ 10"^ 10^ 

m 

Figure 5.8: Solutions of the dispersion relation in the presence of betatron oscil- 
lations in the limit krVth <^ 1, 7 = 30, C = 0.02. Points which do not satisfy the 
inequalities m ^ 1, m^^^Vth < 1 and |Aa;p < Fgy^f^ are plotted in gray. 




Figure 5.9: Fq as a function of m for various values of Vth and the squared ratio of 
the growth rates from Fig. 15.31 and Fig. 15.81 (dashed hne) 
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where ^ = r/rQ and the integration is over the thickness of the layer. The Bessel 
function can be expressed approximately in term of an Airy function as done before. 
We take the linear approximation to the Airy function as discussed previously, and 
this gives 

4/3 



(5.100) 



where C2 ~ 0.259. This is valid for sufficiently big values of 7 and low m. The 
largest values occur for krVth <^ 1, where this quantity is simply 4t>5j. This is 
enough motivation for us to work in this limit. Thus, 

/ 2 \ 

< 27rLiuriclvl \5J^of - ■ (5.101) 

Because we calculated our growth rates in the thin approximation for k ^ it is 
consistent to use 50 = in'^VthV^^ ZrQ6J(j,o. Furthermore, we set u m(j). This is 
consistent even for large growth rates since the exponential growth has stopped. 
With our expression for the saturation amplitude we obtain 

P^<^^jL(lX"'±(^y . (5.102) 

Svr^roe^ |Zp \m J \ 4' / 

Since the number of particles is proportional to ( and the growth rates are 
proportional to v^C fo^' > ^1 the radiated power scales like A^^. This suggests 
that the emitted radiation is coherent. In Fig. 15.101 we plotted the radiated power 
in arbitrary units having evaluated Fq numerically. For large m the curve scales as 
^-5/3 Analytically we obtain with our second approximation for Fq the scaling 
m"3(m^/3/m^/4)4 = m~'^^^. With |Zp ^ 4cj (2/m)^^^ we obtain 



< 3.71 X lO'^m-'- f ^ . (5.103) 
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Figure 5.10: Radiated power (m~^(a;j/0)^) for 7 = 30 and C = 0.02 in arbitrary 
units. The straight hne is proportional to m~^/^ and is shown for comparison. 
Points which do not satisfy the inequahties m » 1, m^^^Vth < 1 and | Acup < Fqv"^^ 
are plotted in gray. 
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5.10 Brightness Temperatures 

We consider the brightness temperatures for conditions relevant to the radio 
emissions of pulsars. Using the Ray leigh- Jeans formula B^, — 2kBTB{i' / cf for the 



where fc^ is Boltzmann's constant and A = 27rrQL is the area of the E- layer. 
The solid angle of the source seen by a distant observer has been computed in 
appendix A and its value is AO = 47r^ro/ (mL). It is assumed that the angular 
size of the source is small such that that radiation from the top and the bottom 
emitted at an angle 9 with respect to the normal is received by the observer at 
the same position. For the sample values 7 = 1000, ( = 0.08, Vth = 0.047~^, 
L — 100 km and m — mi our model predicts a maximum brightness temperature 
of Tb 2 X lO^^K. According to our results from previous sections there may be 
degeneracy from modes with non-zero axial wavenumbers kz- It is reasonable to 
assume that this will increase the brightness temperature by a factor in the order 
of m tan ipcr- Beaming along the z-axis may increase the brightness temperature 
and the observed frequency even further. 

5.11 Applications in Accelerator Physics 

The next-generation linear collider requires a beam with very short bunches and 
low emittance. That is, the beam must occupy a very small volume in phase space. 
The emittance of the pre-accelerated beam is reduced in a damping ring which is 
operated with longer bunches to avoid certain instabilities. The bunch length has 



radiated power per unit area per sterradian at a frequency u = m0/27r gives 




(5.104) 
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to be decreased in a so-called bunch compressor before the beam can be injected 
into the linear collider. A bunch compressor consists of an accelerating part and 
an arc section. Since the bunch lengths of the proposed linear colliders are in the 
order of the wavelength of the synchrotron radiation which is being radiated in the 
arc section, instabilities due to coherent synchrotron have to be taken seriously. 
For a design energy of 2 GeV and 7 x 10^^ electrons per 100 /im our dimensionless 
quantities become 7 = 4000 and ( = 0.08 [70]. Our qualitative analysis of the 
betatron motion suggests that GSR is suppressed for a minimum energy spread of 
vth > C"'7"' = 12.57-^ 

5.12 Discussion and Conclusions 

This work has studied the stability of a collisionless, relativistic, finite-strength, 
cylindrical electron (or positron) layer by solving the Vlasov and Maxwell equa- 
tions. This system is of interest to understanding the high brightness temperature 
coherent synchrotron radio emission of pulsars and the coherent synchrotron ra- 
diation observed in particle accelerators. The considered equilibrium layers have 
a finite 'temperature' and therefore a finite radial thickness. The electrons are 
considered to move either almost perpendicular to a uniform external magnetic 
field or almost parallel to an external toroidal magnetic field. A short wavelength 
instability is found which causes an exponential growth an initial perturbation of 
the charge and current densities. The periodicity of these enhancements can lead 
to coherent emission of synchrotron radiation. Neglecting betatron oscillations 
we obtain an expression for the growth rate which is similar to the one found 
by Goldreich and Keeley [21] if the thermal energy spread is sufficiently small. 
The growth rate increases monotonically approximately as m^^^, where m is the 
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azimuthal mode number which is proportional to the frequency of the radiation. 
With the radial betatron oscillations included, the growth rate varies as m}/^ over 
a significant range before it begins to decrease. 

We argue that the growth of the unstable perturbation saturates when the 
trapping frequency of electrons in the wave becomes comparable to the growth 
rate. Owing to this saturation we can predict the radiation spectrum for a given 
set of parameters. For the realistic case including radial betatron oscillations we 
find a radiation spectrum proportional to m~^/^. This result is in rough agreement 
with observations of radio pulsars [SUIT] (Fig. 13.21) . The power is also proportional 
to the square of the number of particles which indicates that the radiation is 
coherent. Numerical simulations of electron rings based on the fully relativistic, 
electromagnetic particle-in-cell code OOPIC [Ylj recovers the main scalings found 
here. 



Chapter 6 

Particle in Cell Simulations^ 

6.1 Introduction 

Attempts to extract the nonlinear evolution of a plasma are usually unsuccess- 
ful except in some very special cases and numerical methods have to be applied. 
Seeking a straightforward numerical solution of the Vlasov equation however is 
prohibitive except in lower dimensional models |4j because the dimensionality of 
the problem is doubled in a framework which makes use of phase space. A pos- 
sible solution to this problem is MHD which condenses the full momentum space 
distribution to only a few macroscopic quantities like density and current. Nu- 
merical MHD is extremely popular and a vast amount of literature exists on this 
topic. Despite its popularity MHD has some shortcomings. First, the results are 
only as good as the used closing condition. Second, some important effects like 
Landau damping rely on the knowledge of the momentum distribution. Landau 
damping [72] is a stabilization mechanism which is crucial for the generation of 
stable beams in particle accelerators. Therefore, MHD is of limited use in particle 
accelerator physics. Particle tracking programs avoid both problems. The grid 
can be set up in position space and position and momentum can be stored for 
each particle. On contemporary computers this is efficient even for a large number 
of particles. Not working in the limit — > oo anymore the effect of changing 
the number of particles being tracked needs to be investigated. Even though the 
number of particles in a plasma is finite it is usually not possible to track that 

"''This chapter will appear as a journal article |7Tj. Reprinted in modified form 
with kind permission from the American Physical Society, (c) 2005 by the Amer- 
ican Physical Society 
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many and one resorts to tracking N^ac "macroparticles" with each macroparticle 
representing N/Nmac particles. One can only hope that the results obtained with 
Nmac macroparticles are sufficiently close to convergence, i.e. N — > oo, giving a 
reasonably good estimate of the behavior of the real system. 

The objective in this chapter is to simulate the evolution of a particle distribu- 
tion which resembles configuration a in chapter [51 

6.2 Particle-in-Cell Simulations with OOPIC 

For the simulation the software package OOPIC [73] was used. OOPIC is a rel- 
ativistic two-dimensional particle-in-cell code which supports both plain {x, y)- 
geometries and cylindrical (r, z)-geometries. Since the interesting dynamics takes 
places in the azimuthal direction one can only simulate a thin ring (instead of 
a cylinder) in the (x, y)-mode. Loading the initial circular particle distribution 
in the (x,?/)-mode required modifying the source code (files load.cpp, diagn.cpp 
and c_utils.c) to allow the program to handle circular particle distributions. Some 
minor modifications were necessary in order to compile XOOPIC-2.5.1 with gcc 
3.2.2 and the compiler compiler bison 1.28 under SunOS 5.9. The built-in function 
parser was extended to support elliptic integrals. 

Since a thin ring of particles is simulated instead of a cylinder thereof all fields 
and charges were divided by the length L of the cylinder whereas the electron 
mass needs to be divided by L^. L = lOro is chosen unless noted otherwise. The 
electric and magnetic self-fields for a thin ring equilibrium differ from what was 
used in the model. The fields can be found in [71]. It is ensured that OOPIC uses 
these self-fields before the perturbation starts to build up. As it turns out choosing 
the correct self-fields is not too crucial. Leaving them out the system will build 
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them up itself. Once the self-fields are created the system shows no difference in 
behavior. The absence of the self-fields in the dispersion relation might help to 
understand this feature. As in [51] a Gaussian number density profile with RMS 
width Vth was chosen for the initial distribution. 5000 macro particles were tracked 
on a grid with resolution 512 x 512 unless noted otherwise. Once an energy for a 
particle has been chosen it is placed at the equilibrium radius tq = m'-fc{eB)^^, i.e. 
neglecting betatron oscillations particles on the same orbit have the same energy. 
This fixes the azimuthal component of the canonical angular momentum. The 
system can pick up transverse motion quickly. The grid represents a rectangular 
region 40m x 40m big where the ring with radius tq = 10m is centered. 

In Fig. 16.11 the initial particle distribution (gray) and the particle distribution 
after 23ns are shown. The parameters are ( = 0.010, 7 = 30, and Vth = 0.002. 
Qualitatively, a bunching of the particle distribution can be observed. An enlarge- 
ment of a small section of Fig. 16. H is also shown in the same figure. In Fig. 16.21 the 
bunching is shown for successively higher energy spreads. With increasing energy 
spread the bunches become fuzzier and the clean gaps between bunches that can 
be observed for small energy spreads are populated with "stray particles". This 
suggests that it may be harder to achieve complete coherence for larger values of 
Vth- The decoherence due to the non-zero width of the particle beam is investigated 
quantitatively later in the paper. These qualitative features are independent of 7. 

Also note that during the evolution of the circular charge distribution both the 
radius and the width of the ring increase slightly. The former is due to a.) particles 
losing energy and b.) the perturbed magnetic field changing significantly. It tends 
to decrease for small energy spreads and increase for larger energy spreads. Starting 
with a larger radius the radius increases even further, i.e. this is not a relaxation 
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from a "false" to a true equilibrium. Since the non-zero mesh size imposes an 
upper limit on the azimuthal mode number m which can become unstable, it is 
expected that the distance between bunches decreases as the resolution increases. 
This is indeed the case. For larger energy spreads the bunching of the distribution 
becomes hardly visible, but it still can be observed in the z-component of the 
magnetic field (Fig. 16.51) . 

The bunches are slightly tilted and may be connected by a very thin inner ring 
of particles for sufficiently high beam currents. For these reasons it is not possible 
to Fourier transform the charge perturbations in order to compute the growth rates 
for each value of m. Since the resolutions used were low the range of m values is 
restricted. Therefore, only the radiated power is computed which can be obtained 
easily. 

Estimates of the growth rate are two orders of magnitude higher than what 
would be expected. A possible explanation is that the ratio between the saturation 
amplitude and the electric self field 



is typically in the order of 10~^ for the given sample cases which is rather small. 
The initial perturbations due to discreteness, numerical noise etc. are usually in 
the same order of magnitude. Therefore, one cannot expect to see the regime 
covered by the linearized Vlasov equation. This is another reason for focusing 
entirely on the emitted power. 



SEsat 



) 



2 



(6.1) 
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Figure 6.1: Initial particle distribution (gray) and the same distribution after 23ns 
have elapsed. Parameters: ( — 0.010, 7 = 30 and Vth — 0.002 
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Figure 6.2: Particle distribution (7 = 30, C = 0.01) after 23ns for vth = 0.002, 
Vth — 0.008, Vth — 0.015 and Vth — 0.033 (from left to right, top to bottom). All 
lengths are measured in units of meters. 
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Figure 6.3: Particle distribution (7 = 30, Vth = 0.002) after 23ns for ( = 0.005, 
C = 0.010, C = 0.020 and C = 0.040 (from left to right, top to bottom). All lengths 
are measured in units of meters. 
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Figure 6.4: Particle distribution (( = 0.01, Vth = 0.002) after 23ns for 7 = 10, 
7 = 30, 7 = 75 and 7 = 90 (from left to right, top to bottom). All lengths are 
measured in units of meters. 




Figure 6.5: z-component of the magnetic field (self-field plus perturbation without 
external magnetic field) after 23ns for ( — 0.010, 7 = 30 and Vth — 0.025. The size 
of the area depicted is 40m x 40m. 
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Figure 6.6: Loss of kinetic energy in W vs. 7 for ^ = 0.01 and Vth = 0.002. 
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Figure 6.7: Loss of kinetic energy in W vs. ( for 'j — 30 and Vth — 0.025. The 
solid line shows the best fit. 
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6.3 Radiated Power 

In Fig. 16.61 and Fig. 16.71 the radiated power determined by measuring the kinetic 
energy loss of the electron cloud after approximately 2.36 ns is plotted as a function 
of 7 and (, respectively. After 0.24 ns the perturbations have saturated and the 
emitted power is fairly constant. A quadratic dependence can be established, i.e. 
the first two relevant scalings expected from the analytical model are recovered. 

The simulation has been repeated at lower (256 x 256) and higher (1024 x 
1024) resolution. No significant effect could be observed. This is consistent with 
the model which predicts that most power is emitted by modes with low m. The 
resolution is not high enough to resolve the path length difference of orbits with 
different radii at these low values of Vth- Also, decreasing the stepsize dt to 2.5 ps 
and increasing the number of macro particles to 50000 has a negligible effect. 

Finally, the effect of the energy spread is investigated. With increasing Vth the 
power decreases which is due to the decoherence described by the factor Fq defined 
in Eq. ( 15. 88^ . The results are plotted in Fig. l6.8l for the parameters ( = 0.01, 7 = 30 
and ( = 0.005, 7 = 10, respectively, and L = vq. In the former case mi is 27 and 
in the latter case it is 0.4. Eq. (15.1031) becomes 

P^3.71 X 10^^^L^^y^m[t7rJm{ujro)Hm{ioro)f (6.2) 

Despite mi being much larger than 1 for the first set of parameters the slopes 
in Fig. 16.81 match exactly only if the summation starts at m = 1. This suggests 
that modes with m < mi do radiate and can be described by the same disper- 
sion relation. Since the power scales as m~^/^ these modes may actually be very 
important for computing the total energy loss. Note that while the simulation 
suggests P oc Eq. (I5.103P (which was derived under the assumption L > r^) 
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gives P (X L. A 2D simulation cannot explain how the radiation from different 
axial positions on the cylinder interacts. In the thin ring case doubling L doubles 
the number of particles and therefore quadruples P. Fortunately, as can be 
seen in the derivation of Eq. fl5.103p in [51] the C, Vth and 7 dependent part of P is 
independent of L. In Fig. 16.81 the overall factor matches if tq = lOOL whereas the 
growth rate for a perturbation of a cylinder and a thin ring coincide for tq = L 
[5T] . Also note that Fig. 16.11 suggests kr-Vth^o ~ 1, whereas Eq. (15.98^ was derived 
under the assumption krVthTo -C 1. 

6.4 Conclusions 

The particle in cell code OOPIC was used to simulate the evolution of density 
perturbations in a thin ring of charged particles which move in relativistic almost 
circular motion in an external magnetic field. The results were compared with the 
model in [51]. Comparisons of the simulation with the model shows approximate 
agreement with the main predicted scaling relations. In particular the bunching 
effect could be observed very clearly and the emitted power is proportional to the 
square on the number density which implies coherent radiation. The dependence 
on the energy spread can be recovered exactly assuming all modes contribute to 
the observed energy loss suggesting that the model may apply even if m < mi. 
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Figure 6.8: Total power radiated as obtained from OOPIC for the parameters 
( = 0.01, 7 = 30 (solid) and C = 0.005, 7 = 10 (dashed). The dash-dotted hne is 
proportional to Xlm=i ^rn- 



Chapter 7 

MHD Approach for a Brillouin Flow 

7.1 Theory 

We consider a laminar Brillouin type equilibrium of a long, non-neutral, cylindrical 
relativistic electron (or positron) layer in a uniform external magnetic field Be = 
BgZ, where we use a non-rotating cylindrical (r, 0, z) coordinate system. The 
electron velocity is v = v^{r)4> = v{r)4> = r(j){r)4> The self-magnetic field is in 
the 2;— direction while the self-electric field is in the r— direction. The radial force 
balance of the equilibrium is 

-70V = — + , (7.1) 
me 

where 7 = (1 — t>^)~^/^ is the Lorentz factor with velocities measured in units of 
the speed of light, B = 3^, + Bgis the total (self plus external) axial magnetic field, 
E is the total (= self) radial electric field, and q and me are the particle charge 
and rest mass. We have 

ld{rE) ^ dB ^ 

— -r^ = 47rpe , -r = -^TiPeV , (7.2) 
r ar ar 

where Pe{r) is the charge density of the electron layer. 

We consider weak layers in the sense that the 'field reversal' parameter 



r2 



C = - — / drpeV (7.3) 



ri 



Be 

is small compared with unity, ("^ <^ 1. Under this condition Eq. (17. ip gives = 
—qBf./{me^). Here, we have assumed that the layer exists between ri and r2. We 
also consider that the Lorentz factor is appreciably larger than unity in the sense 
that 7^ > 1. 
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We consider general electromagnetic perturbations of the electron layer with 
the perturbations proportional to 



/a(r) exp{im(f) — iut) 



(7.4) 



where a = 1, 2, .. for the different scalar quantities, m = integer, and u the angular 
frequency of the perturbation. Thus the perturbations give rise to field components 
6Er, SEfj,, and SB^. The perturbed equation of motion is 

d 



dt 



(v + 5v) ■ V 



(7V + v57 + 7(5v) 



= — (5E + V X 5B + 5v X B) , (7.5) 

me 

where the deltas indicate perturbation quantities. This equation can be simplified 
to give 



-i-fAto -70(1 + 7^) - ^5 

7(^ + (7(^r)' + -i-f^Au 



SVr 



A. 

me 



5Er + v6B^ 



(7.6) 



where the prime denotes a derivative with respect to r, and 

Aco'(r) = uj — m(f){r) 

is the Doppler shifted frequency seen by a particle rotating at 0, 

Using the equilibrium equation (17. ip and the condition <^ 1, the matrix in 
Eq. (17.61) is approximately 



V 



—i'yAuj —7^0 
7^(0r)' —i'y^Auj 



(7.7) 
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We have used the fact that {'^(pr)' = 7'^(0r)'. For ^ 1 we have (0r)' = 0/7^ 
and 4>' = —v'^4>/r. Consequently 



det(P) = 7^(02 - Auj^' 



Inverting Eq. (I7.6p gives 



SVr 



and 



medet (V) 



QT 
medet(P) 



iAuj{5Er + v5Bz) + 



{7.t 



(7.9) 



(7.10) 



7.2 Field Sources 



The source terms due to the perturbation are 



SJr = Pe^Vr , ^ud 6J^ = PeSv^ + Sp^V 



and from the continuity equation, 



1 



Spe = [DriPeSVr) + ik^PeSv^] , 



(7.11) 



(7.12) 



where Dr = {l/r)[d/dr{r...)] and = m/r is the azimuthal wavenumber. Ex- 
panding this equation gives 

1 r dj" 



Spe 



iAuj 



+ J'[- ik^{5Er + v5B,) + a,{5E^ 



Here, 



QPe 



medet(X') 

has the role of the distribution function of angular momentum [75] . 



(7.13) 



(7.14) 
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7.3 The Limit Au <C 7 ^ 



For 



< 1 



the resonant term in 5pe proportional to l/Au is dominant. 
Thus we have 



Pe^V^, = — -{6Er + v5B,) + C ( 



r 



6peV 



rcf) 
iAu 



dr 



6E. + .. 



+ 



Aoj 



(7.15) 



where the eUipsis indicates a term equal to the middle line of equation fl7.13l) . To 
leading order in Auj\ we have 

5Jr = , 



or 



SJa, 



rcj) 
iAu 



dr 



6E^ + J^[- ik^{5Er + v5B,) + Dr{5E^)\ 



In this approximation we also have 



Pe 

B 



For an electron layer with i?e > 0, we have (p > and J-' > 0. 
Thus, 



iq 



Sv^ -- 

where ACo = Auo / {m(p) . Finally, 



—niAuj 



m'j'^Auj 



SE, 



iq 

nie'^^Auj 



6E. 



(7.16) 



(7.17) 



(7.18) 



(7.19) 
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The reason for the discrepancy from Eq. ( I5.40p is the difference in the used equi- 
hbrium. In chapter [5] the non-zero width was caused by betatron oscillations of 
particles with the same average angular velocity. In a Brillouin flow particles on 
different orbits have different angular velocities. This additional source of shear 
is reflected in the (70r)' term. Setting this term to zero one obtains the previous 
results from chapter [5] again (cf. Eq. (17.221) ). 

In the absence of radial currents the linearized continuity equation simply reads 

5p = -l^PoSv^ (7.20) 



Thus, 



Sp= 'f^jxh^^^ ^^-21) 



7.4 The Limit Aw > 

If we completly neglected radial motion as it was done in an earlier paragraph in 
chapter O Eq. fl7. 101) would read 

H = —^5E^ . (7.22) 

In our two-dimensional MHD approach the motion is not constrained to a 
fixed radius and we are wondering under which conditions the radial motion in 
an unconstrained model can be neglected, i.e. when Eq. (I7.22p and Eq. (17.101) 
coincide. 5vr = implies that the forces due to 6E^, 6Bz and 6Er have to balance. 
Eq. dH]) gives 

(j)5E^ = iAuj{5Er + v5B,) (7.23) 
Thus, Eq. (03) and Eq. fTTTUD coincide if 

Acu > (7.24) 
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or if we choose an equilibrium distribution with zero average shear. In both cases 
the growth rates are given by Eq. (15.441) . 



7.5 Configuration b 

In this section we are going to investigate the stability properties of an equilibrium 
with the same number density and velocity profile as before, but with different 
external fields. Instead of an external magnetic field in the z direction we consider 
an equilibrium with an azimuthal magnetic field acting as a guiding field and 
a radial electric field. The latter is included in the equilibrium condition and 
therefore does not enter the linearized Euler equation. would only enter if we 
considered motion in the axial direction and non-zero axial wavenumbers. Thus, 
we obtain the matrix T> again without the Bq terms, i.e. for 7 ^ 1 



V 



-i'-fAuj —7^0 
270 —i'y^Aui 



(7.25) 



with 



We obtain 



det(P) =7^(20^- Acu') 



6Vr 



and 



medet (V) 



- iAuj{6Er + v6B^) + 



(7.26) 



(7.27) 



(7.28) 



medet(P) 

In the limit {Auj)"^ ^ 20^ the azimuthal current is given by Eq. (17.221) again and 
by Eq. (mm for (Acu)^ < 20^. 
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7.6 Two Cylinder Model 

Instead of solving the full two-dimensional problem we solve it for two concentric 
cylinders, i.e. for the number density we have 

n{r) = ^noV2Trvthro {S{r - ri) + 6{r - r2)) (7.29) 

Some interesting results can be obtained in that limit assuming zero average shear. 

The linearized continuity equation becomes 

1 d %Tfi 
- icuSp + {rpoSvr + rvorSp) + — {poSv^ + Vo^Sp) = (7.30) 

We drop all radial derivatives in the continuity equation which is consistent with 
our previous approximation. Thus, 

Sp = '-e'noV2^Vtnrl {S{r - n) + S{r - r2)) (7.31) 
The Green function is 

5$(ri) = 27r^i {riJm{uJn)H^^{u;n)5pi + r2Jm{^r2)H^\u;n)6p2) (7.32) 
5$(r2) = 27r^t {nJ^{ujr2)Hg\ujn)6p^ + J^(a;r2)ifi^)(^^2)5p2) , (7.33) 

so the problem can be written as a matrix 



1 1 

A = -27r^ie'^-noV27rvthrl—'y~'^ 
Z rl 



(7.34) 



y (Aa;(ri)) ^ii J^(a;ri)i/„(cjr2) (Acj(r2)) '^^/m{ujr2)H^{ujr2) J 
acting on the column vector ((5$(ri)), S^{r2))^ which returns the same column 
vector. Nontrivial solutions can only exist if 

det(A - 1) = . (7.35) 
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Finally, this equation can be solved for Acj. 

The coefficient in front of the matrix is equal to — |iC7~^- We introduce the 
notation Au; = Acj(ro), Acj(r) = (1 + Acj)^ — 1 and solve the dispersion relation 
for Ad; approximating the Bessel functions by Airy functions 

Jm{z) ={-] Ai {-{2/mf\z - m)) (7.36) 

Y^{z) = -{-] Bi {-{2/mf/\z - m)) (7.37) 

and neglecting terms proportional to AcD which are assumed to be small compared 
with 7~^. The growth rates for the parameters C, = 0.02, 7 = 30, m = 500 are 
plotted in Fig. 17. 1[ 

Numerically, we find that the drop occurs very roughly for Vth = jir/m if j is an 
even integer, i.e. when the thickness of the layer is an integer multiple of 4 times 
the wavelength of the radiation. This periodicity is due to the Bessel functions. 
The drop in the growth rate is caused by the Coulomb term as can be seen by 
using the following approximation for the Airy functions 

Ai{w) ^ ci - C2W (7.38) 
+ C2w] (7.39) 

where Ci = l/[32/3r(2/3)] and C2 = l/[3^/^T{l/3)] which is justified for jwp < 
0.5. This gives Z^(w) = iJ^{w)H^{w) ^ {2/m)'^/^[VS{cl - 4w) + i{ci - caw)^]. 
Evaluating the Green function at ro for w = we just obtain the radiation term 
quoted by Goldreich and Keeley [21] and the two unstable modes shown in Fig. 17.11 
become degenerate. Keeping the ffist order term in w we recover the Coulomb term 
as well. This approximation is valid in the shaded area of Fig. I7.1[ 

The dependence of the larger mode in Fig. 17.11 on Vth seems to be in rough 
agreement with the dependence obtained from including the decoherence due to 
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Figure 7.1: Growth rates as a function of energy spread. The Airy fuctions were 
retained. Parameters: C = 0.02, 7 = 30, m = 500 
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betatron oscillations in chapter [5l We speculate that the betatron oscillations 
themselves are not too important. What is important is the fact that particles 
move on different orbits with different angular velocities 0(r). 



Chapter 8 
Summary 

What have we learnt in the last three chapters? We analyzed the stability prop- 
erties and the power spectra of charged particles executing circular motion at rel- 
ativistic speeds using three different techniques (Vlasov, PIC, MHD). Relativistic 
plasmas are capable of self-bunching and emitting coherent synchrotron radiation. 
This has been established by all three methods. The particular field geometry of 
the external guiding fields played only a minor role (MHD). In particular there 
were only minor differences in the stability properties of configuration a (Fig. 15.11) 
and configuration b (Fig. 15. 2p . The self- fields do not enter the dispersion relation 
either as long as C ^ 1 (Vlasov, PIC). For the Vlasov treatment we selected a 
rather special distribution function with no spread in the azimuthal component of 
the canonical angular momentum and a small energy spread. Other choices would 
have been possible. However, since we are able to recover the main effects with 
MHD the only important kinetic effect seems to be the Landau damping observed 
for non-zero axial wavenumbers (Vlasov). We put a lot of effort in understand- 
ing the effect of a small energy spread. In the absence of such an energy spread 
the results by Goldreich and Keeley were recovered [21]. Once the energy spread 
exceeds a critical threshold (which is a function of the azimuthal mode number 
m) the growth rate obtained by Goldreich and Keeley is attenuated by a factor 
which is due to the decoherence caused by the betatron oscillations. Once this 
threshold has been exceeded there is no further dependence of this factor on the 
energy spread as long as the spread remains small. Fig. 15.91 summarizes these 
findings. The characteristic power spectrum which scales as z/~^/^ is due to this 
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decoherence and it is in good agreement with observations from actual radio pul- 
sars. In a Brillouin flow the CSR instability may be almost completely suppressed 
by the shear intrinsic to such an equilibrium. The amount of shear is very small 
and the resolution of the PIC simulations were not high enough to resolve the path 
length differences of orbits with different radii at the low energy spreads that were 
used. In the Vlasov approach we used a different equilibrium where the non-zero 
thickness was due to betatron oscillations. At least on average all particles were 
moving at the same angular velocity. Still, both the naive two cylinder model 
and the Vlasov approach for a thin layer with betatron oscillations give the same 
dependence of the growth rate on the energy spread, but further investigation is 
needed to establish this result. Whether the negative mass instability is present in 
the system under investigation is a tough call. As Fig. 17.11 and the accompanying 
paragraph suggest the Coulomb term tends to stabilize the second unstable mode. 
This conclusion is in agreement with earlier findings by Goldreich and Keeley [21] . 
However, the very presence of a second unstable mode might be an artifact of 
the negative mass instability. The treatment in ^21j and chapter [5] differ in the 
following way: In the former an increase in energy leads to higher angular velocity 
whereas in the latter the angular velocity decreases because of the negative effec- 
tive mass. Still the results resemble each other. In some sense the CSR instability 
in chapter O is a "negative mass instability" with a dominating radiation term 
instead of a Coulomb term. The effect of the Coulomb term can only be seen for 
non-zero energy spreads. In Fig. 17.11 the second mode does drop down to zero for 
larger energy spreads which is a characteristic of the negative mass instability, but 
a slight increase in energy spread leads to a recovery of that mode again. Both 
the negative mass instability and the CSR instability manifest themselves as a 
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bunching of the distribution. It is conceivable that the two instabihties interfere 
with each other in the sense that the weak negative mass instabihty caused by the 
Coulomb term can only disturb the CSR instability caused by the much stronger 
radiation term. In this case it would be very hard to tell them apart. 

As was pointed out in chapter [3] there are a couple of competing theories in the 
literature trying to explain the radio emission of pulsars. With the exception of 
the CSR instability none of them seems viable. The biggest obstacle to applying 
the CSR instability to pulsars has been a lack of more detailed models. We hope 
that we bridged that gap. It is certainly far from easy to figure out how a system 
evolves by only looking at its saturated state. However, our model is not arbitrary. 
It is well motivated and rather simple. Only very few dimensionless parameters 
are needed and its results are strikingly generic. We certainly do not claim that 
our model is the final word, but dismissing it as a coincidence without further 
studies would be too easy. The fact that many of the presented results only 
depend on dimensionless parameters implies scale invariance. This opens up the 
possibility of testing pulsar radiation mechanisms in the lab. Particle accelerators 
suitable for such experiments already exist and the similarity of instabilities found 
in accelerators and astrophysical objects is hard to overlook. Proof of principle 
experiments for different astrophysical problems have already been performed |76j . 
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Green's Function 

The Green's function for the potentials give 
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where G is the Fourier transform of the Green's function. The "C" on the integral 
indicates an cu— integration parallel to but above the real axis, Im(u;) > 0, so as to 
give the retarded Green's function. 

Because of the assumed dependences of Eq. fl5.26l) . we have for the electric 
potential, 
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- {u^ - kl) 

where = fc^ + A;^. Because uo has a positive imaginary part, this solution corre- 
sponds to the retarded field. Also because Im(ci;) > 0, the /t— integration can be 
done by a contour integration as discussed in [77] which gives 
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where k = (u;^ — klY^"^, where r< (r>) is the lesser (greater) of {r,r'), and where 
Hm\x) = Jm{x) +iYm{x) is the Hankel function of the first kind. From the Lorentz 
gauge condition 



ujrnk. 
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Eqs. ( IA.4P and flA.SP are useful in subsequent calculations. 

To determine the total synchrotron radiation from the E-layer it is sufficient to 
calculate 5 A at a large distance from the E-layer. We assume that the E-layer has 
a finite axial length and exists between —L/2 < z < L/2. Thus we evaluate 5 A 
in a spherical coordinate system R = (i?, 6*, 0) at a distance L. The retarded 
solution is 
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im(j)' + ikzz' — iuj{t + H ■ r') 

(see, e.g. ch. 9 of [78j). The source point is at (x' = r'cos0', y' = r'sin0', z'). The 
observation point is taken to be at (x = 0, y = Rsm9, z = Rcos9). Consequently, 
R ■ r' = r' sin 6' sin 0' + z'cos6'. The phase factor exp{iuR) does not affect the 
radiated power and is henceforth dropped. 

For the cases where 6J^ is the dominant component of the current-density 



perturbation we have 



m 

R 



r'dr'c 



sm (f) , cos 



6Jif,{r') exp{im(f)' — iur' sin 9 sin (p') , 



(A.7) 



where 



S{9) = L 



sin[(A;2 — uj cos9)L/2] 
{k,-ujcos9)L/2 



(A. 



128 



is a structure function accounting for the finite axial length of the E-layer, and u 
superscript indicates u = m(f). Carrying out the 0' integration in Eq. ( 1A.7I) gives 

s{e) 
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r'dr' 6J4r')[..] 
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where 
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and where the prime on the Bessel function indicates its derivative with respect to 
its argument. The radiated power per unit solid angle is 
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where k = cijR is the far field wavevector. 

For a radially thin E-layer, (Ar/ro)^ <^ 1, Eqs. (1A.9I) and 



.101) give 
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The factor within the curly brackets is the same as that for the radiation pattern 

of a single charged particle (see ch. 9 of [78]). 

The factor S'^{9) in Eq. (lA.lip tightly constrains the radiation to be in the 

direction 9^, = cos^^ (kz/uj) if the angular width of S'^{9), the half-power half-width 

A9i/2 ~ tt/{ujL), is small compared with the angular spread of the single particle 

synchrotron radiation, I/7, which is the angular width due to the Bessel function 

terms in Eq. (lA.llI) . This corresponds to E-layers with L ^ TC'y /uj = 7rro7/m. For 

L ~ To, we need m ^ nj, which is satisfied by the spectra discussed later in §5.81 

In this case, Eq. (lA.lip can be integrated over the solid angle to give 
71 L sin 9^: 
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One limit of interest of Eq. (1A.12|) is that where fc^ = so that 6^ = 7i/2 and 

2 
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where we have set u mcf). The total radiated power is P = J2m ^rn- 



Appendix B 

Bessel Function Approximations 

The Bessel functions Jm{ojr) and Ymiojr) are the two hnear independent solutions 
of the differential equation 



-^r^- — +A5E,{r) = Q (B.l) 
r or or r J 

Computing these functions is very cumbersome for high values of m - even on 
modern computers. If analytical results are desired approximating those functions 
may be a necessity. The following approximations are extremely useful |79[ [77] 

/ 2 \ 

Jm (m + zm^''^) = f — ] Ai + 0{m-^) 

/ 2 \ 

y™(m + W/3) = - f — j Bi (-21/3^) + O(m-i) (B.2) 

The Airy functions can be evaluated very quickly on modern computers. For 
analytical results further approximations may be required. 

For \w\^ > 1 

ki{w) ^ (2v/^)-iw;-i/^exp(-2^3/2/3) 

Bi(^) ^ (v/^)-iw;-i/^exp(2ti;3/V3) (B.3) 

and for |wp < 0.5 

Ai{w) = ci - C2W + 0{w^) 
Bi{w) = Vs[ci + C2W + 0{w^)] (B.4) 

where ci = l/[32/3r(2/3)] ^ 0.355 and cs = l/[3^/^T{l/3)] ^ 0.259. 
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Appendix C 

Source Code Listings (Maple 
Worksheets) 

The following worksheets were generated using Maple 6 for IRIX 6.5. 
C.l Solver for Eq. 

> eqn : =-H-Pi*zeta*Z* (Doint-l/gainma0"2) /Doint"2 ; 



> Z : = (2/ni) " (2/3) * (AiryAi (w) *AiryBi (w) +I*AiryAi (w) "2) ; 



> w: = (ni/2)~(2/3)*(l/gaiimiaO"2-2*Donit) ; 



> gainma0:=30: zeta:=0.02: 



> m:=1000; 



> f solve(eqn,Doint=le-6+le-6*I) ; 

-0.00001518689566 + 0.0005418589616 I 

> Im(%)*m; 

0.5418589616 
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C.2 Solver for Eq. (jSOT) 

> gainma0:=30: zeta:=0.02: vth:=l/gainmaO: 
psi : =arctan(0 . 005) : in:=le3: eps:=le-4: 

> Z : = (2/ni) " (2/3) * (AiryAi (w) *AiryBi (w) +I*AiryAi (w) "2) : 

> w: = (ni/2)~(2/3)*(l/gainma0"2+tan(psi)"2-2*u*tan(psi)) : 

> F : =unapply (Heaviside (-lm(z) ) *I*sqrt (2*Pi) *exp (-z"2/2) +l/sqrt (2*Pi) 

*Int(exp(-x~2/2)/ (x-z) ,x=-inf inity . . infinity, 
digits=5 ,method=_NCrule) ,z) : 

> C:=unapply(Pi/vth~2*zeta*Z/tan(psi)*(u-l/gainmaO"2/tan(psi)) ,u) : 

> B :=unapply (Pi/vth"3*zeta*Z* (u-l/gainma0"2/tan(psi) )* 

(0*l-u/tan(psi)) ,u) : 

> u: = (0. l+0.3*I)/in/tan(psi) ; 

u := 0.02000000000 + 0.06000000000 I 

> evalf ( 1-B (u) *F (u/vth) +C (u) ) ; 

-0.278474072 - 0.9518667283 I 

> for i from to 8 do 

> epsl :=evalf (l-B(u+eps)*F((u+eps)/vth)+C(u+eps)) ; 

> epsO : =evalf (1-B (u) *F (u/vth) +C (u) ) ; 

> un: =u-eps/(epsl-epsO) *epsO : 

> u:=un; printfC'/.g Ig "/.g Ig \n" ,Re(u) ,Im(u) ,Re(epsO) ,Im(epsO)) ; 

> end: 

> Domt :=u*tan(psi)/(l-l/2/gainma0"2) ; 

Domt := -0.00003303461986 + 0.0004692459188 I 

> Im(Domt) *m; 
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0.4692459188 

> Re(Domt)*m; 

-0.03303461986 



C.3 Evaluator for Eq. ( 157961 ) 

> Ff :=unapply(I~n*exp(-I*n*arctan(krb/m) )/2/Pi* 

Int (exp (-I*n*theta-Chi"2/2* (m/sqrt (m"2+krb"2) -sin(theta) ) ~2) , 
theta=-Pi. .Pi) ,n,Clii) ; 

> gainma0:=30: zeta:=0.02: vth:=l/gainmaO: 

> ni:=5e4; krb:=10: 

> 

> F[0] :=evalf (Ff (0,in*vth)) ; 
> 

> Domt : =(2/Pi) ~ (1/4) *sqrt (-zeta*F [0] ) /sqrt (vth* (m~2+gamma0~2*krb~2) ) ; 

> evalf (ni*Im(Domt) ) ; 



0.07543389577 
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C.4 Solver for Eq. ( IHTOTh 

> eqn:=l=Pi*zeta*Z*(-l/gainmaO"2*sum(F[n]/(Domt-n/l)~2,n=-N. .N)); 



> Ff :=unapply(I"n/2/Pi*Int(exp(-I*n*theta-Chi"2/2*sin(theta)"2 + 

Chi~2*sin(theta)-Chi"2/2) ,theta=-Pi. .Pi) ,n,Chi) ; 

> Z : = (2/1) ~ (2/3) * (AiryAi (w) *AiryBi (w) +I*AiryAi (w) "2) ; 

> w: = (l/2)~(2/3)*(l/gaiimia0"2) ; 

> N:=0: 

> gammaO : =4000 : zeta:=0.08: vth:=0.04: 

> l:=le3; 

> for j from -N to N do 

> F[j] :=evalf (Ff (j ,l*vth)) ; 

> od; 

> solve (eqn, Domt) ; 

> l*Im(%[2]); 

0.002124396570 



Appendix D 



Source Code Listings (XOOPIC) 

The results presented in chapter [6] were obtained using a modified version of 
XOOPIC-2.5.1. In the following sections we present the output of UNIX diff 
command applied to the original and the modified source files. Together with the 
original source files it is possible to recover the modified version of XOOPIC. 

D.l File c_utils.c 

> diff /tmp/oopic/otools/c_utils . c ~/oopic/otools/c_utils.c 
61a62 

> /* wrappers for atan2 and erf 
69a71 

> float ATAN2W (double y, double x) i return atan2(y, x) ; } 
70a73,74 

> float ERFW(double x) {return erf (x) ; } 
> 

92a97,115 

> /* some useful functions 
> 

> float EllipticE(float x) { 

> float ml, res; 

> ml = 1-x; 

> res = l+.46301*ml+. 10778*ml*ml + 

( . 24527*ml+ . 04124*ml*ml) *log(l/ml) ; 
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> if (x<0.999) return res; 

> return 1 . ; 

> } 
> 

> float EllipticK (float x) { 

> float ml, res; 

> ml = 1-x; 

> res = 1. 38629+. 11197*ml+.07252*ml*ml + 

(.5+. 12134*ml+.02887*ml*ml)*log(l/ml) ; 

> if (x<0.999) return res; 

> return 4.5; 
> 

> } 
> 

125dl47 
< 

D.2 File dump.cpp 

> diff /tmp/oopic/otools/dump . cpp ~/oopic/otools/dump . cpp 
Oal 

> extern "C++" void write_validation() ; 
58a60,65 

> 
> 

> // 
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> 
> 

> /* 

72a80,81 
> 

> */ 

73a83,86 

> // 

> // 
> 

> write_validation() ; 
75d87 

< 

D.3 File diagn.cpp 

> diff /tmp/oopic/otools/diagn. cpp "/oopic/otools/diagn. cpp 
879c879,881 

< #ifdef BENCHMARK 

> // #ifdef BENCHMARK 

> // 

> // 
885c887 

< void write_validation() { 
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> void write_validation2() { 
888a891,900 

> int j , k ; 

> double Bval [1024] , Btr [1024] ; 

> double Eval [1024] , Etr [1024] ; 

> double rhoval [1024] ,rhotr [1024] ; 

> double a_sum,b_suin; 

> double R,max, count, count mcix,mcixR; 
> 

> 
> 
> 

890,892d901 

< fordnt j=0;j<J;j++) 

< fordnt k=0;k<K;k++) 

< fprintf(trace_f lie, "7.10. 4g\n",E[j] [k].el()); 
893a903,975 

> for (R=110;R<=135;R+=5) 

> { 

> max = 2*3.1415*R; 

> f or(j=0; j<max; j++) 

> { 

> Bval[j] = (theSpace->getBNodeDynamic() ) [(int) 

(J/2+R*sin(j/R))] [(int) (J/2+R*cos(j/R))] .e3() ; 

> fprintf (trace_file,"y.d %d '/.d %10.4g\n" ,0, (int) R, j ,Bval [j] ) ; 
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} 

f or ( j =0 ; j <max ; j ++) 
{ 

a_sum = 0; 
b_sum = 0; 
for (k=0 ; k<max ; k++) 
{ 

a. sum += Bval [k] *cos(j*2*3. 1415*k/max) ; 

b. sum += Bval [k] *sin(j*2*3. 1415*k/max) ; 
} 

a_sum = 2*a_suni/max; 
b_sum = 2*b_sum/max; 

Btr[j] = sqrt (a_sum*a_sum+b_suiii*b_sum) ; 

fprintf (trace_file,"y.d %d %d y.l0.4g\n" , 1 , (int) R,j,Btr[j]); 
} 



for ( j =0 ; j <max ; j++) 
{ 

EvalEj] = (theSpace->getENode()) [(int) ( J/2+R*sin(j/R) )] 
[(int) (J/2+R*cos(j/R))] .el()*(-l)*cos(j/R)+ 

(theSpace->getENode()) [(int) ( J/2+R*sin(j/R) )] 
[(int) (J/2+R*cos(j/R))] .e2()*sin(j/R) ; */ 

Eval[j] = sqrt(pow((theSpace->getENode()) 
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[(int) (J/2+R*sin(j/R))] 

[(int) (J/2+R*cos(j/R))] .el() ,2)+ 

> pow ( (theSpace->getENode ( ) ) 

[(int) (J/2+R*sin(j/R))] 

[(int) (J/2+R*cos(j/R))] .e2() ,2)) ; 

> fprintf (trace_file,"y.d 7A Id %10.4g\n" ,2, (int) R, j ,Eval [j] ) ; 

> } 
> 

> f or(j=0; j<max; j++) 

> { 

> a_sum = 0; 

> b_sum = 0; 

> for(k=0;k<max;k++) 

> { 

> a.sum += Eval [k] *cos ( j *2*3 . 1415*k/max) ; 

> b_sum += Eval [k] *sin(j*2*3. 1415*k/max) ; 

> } 

> a_sum = 2*a_suin/max; 

> b_sum = 2*b_suin/max; 

> Etr[j] = sqrt (a_suin*a_sum+b_suin*b_suin) ; 

> fprintf (trace_file,"y.d "/.d '/A y.l0.4g\n" ,3, (int) R,j,Etr[j]); 

> } 
> 

> f or(j=0; j<max; j++) 

> { 
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> rhovalEj] = (theSpace->getRho() ) 

[(int) (J/2+R*sin(j/R))] [(int) (J/2+R*cos(j/R))] ; 

> fprintf (trace_file,"y.d %d %d %10.4g\n" ,4, (int) R, j .rhoval [j] ) ; 

> } 
> 

> f or(j=0; j<mcLx; 

> { 

> a_siim = 0; 

> b_sum = 0; 

> for(k=0;k<max;k++) 

> { 

> a.sum += rhoval[k]*cos(j*2*3.1415*k/max) ; 

> b.sum += rhoval[k]*sin(j*2*3.1415*k/max) ; 

> } 

> a_sum = 2*a_suin/max; 

> b_sum = 2*b_siain/mcLx; 

> rhotr[j] = sqrt(a_suin*a_sum+b_sum*b_suin) ; 

> fprintf (trace_file,"y.d %d Id y.l0.4g\n" ,5, 

(int) R, j ,rhotr [j] ) ; 

> } 
> 

> } 
> 

> 



896c978,1012 
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< #endif 
> 

> void write_validation() { 

> FILE *trace_file; 

> int J = theSpace->getJ() ; 

> int K = theSpace->getK() ; 

> int j , k ; 

> double minB,maxB,minE,maxE,temp; 
> 

> 

> if ((trace_file=fopen("trace.dat","a"))==NULL) exit(l); 
> 

> 

> minB = (theSpace->getBNodeDynamic() ) [0] [0] .e3() ; 

> maxB = (theSpace->getBNodeDynainic()) [0] [0] .e3() ; 
> 

> minE = sqrt(pow((theSpace->getENode()) [0] [0] .el() ,2) + 

pow((theSpace->getENode()) [0] [0] .e2() ,2)) ; 

> maxE = sqrt(pow((theSpace->getENode()) [0] [0] .elO ,2) + 

pow((theSpace->getENode()) [0] [0] .e2(),2)); 

> 
> 

> for(j=0;j<J;j++) 

> for(k=0;k<K;k++) 
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> { 

> if ((theSpace->getBNodeDynainic()) [j] [k] .e3() < minB) 

minB = (theSpace->getBNodeDynaiiiic()) [j] [k] .e3() ; 

> if ((theSpace->getBNodeDynainic()) [j] [k] .e3() > maxB) 

maxB = (theSpace->getBNodeDynainic()) [j] [k] .eSO ; 

> temp = sqrt(pow((theSpace->getENode()) [j] [k] .el() ,2) + 

pow((theSpace->getENode()) [j] [k] . e2() ,2) ) ; 

> if (temp < minE) minE = temp; 

> if (temp > maxE) maxE = temp; 

> } 

> fprintf (trace_file,"y.l0.4g 7.10. 4g yol0.4g 7.10. 4g\n", 

minB, ma:xB, minE, maxE) ; 

> 

> f close(trace_f ile) ; 
> 

> } 
> 

> // #endif 

D.4 File evaluator.y 

> diff /tmp/oopic/otools/evaluator .y ~/oopic/otools/evaluator .y 
83a84 

> ATAN2W(); 
87a89,90 

> EllipticEO ; 
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> EllipticKO ; 
91a95 

> ERFWO; 
97al02 

> "ataii2",ATAN2W, 
104allO,lll 

> "EllipticE",EllipticE, 

> "EllipticK",EllipticK, 
105all3 

> "erf", ERFW, 

D.5 File evaluator.h 

> diff /tmp/oopic/otools/evaluator .h ~/oopic/otools/evaluator .h 
42a43 

> float ATAN2W(float y, float x) {return (f loat)atan2(y,x) ; } 
47a49 

> // float ATAN2W(float ); 

D.6 File load.cpp 

> diff /tmp/oopic/physics/load. cpp ~/oopic/physics/load. cpp 
146al47 

> Scalar Jinx,Kmx; 
211a213 



145 



213c215,222 

< u = incLxwellian->get_U() ; 

> Jmx=grid->get J() ; 

> Kmx=grid->getK() ; 

> u = maxwellian->get_vO() ; 

> Vectors beta = iSPEED_OF_LIGHT*u; 

> Scalar gammaO =l/sqrt(l-beta*beta)*sqrt((x-Vector2(Jmx/2,Kmx/2)) 

* (x-Vector2 ( Jmx/2 , Kmx/2) ) ) / ( Jmx/4) ; 

> Scalar phi = atan2(x*Vector2(l,0)-Kmx/2,x*Vector2(0, 1)- Jmx/2) ; 

> u = -gainma0*SPEED_0F_LIGHT*(l-l/2/gaimna0/gainma0) 

*cos(phi)*Vector3(l,0,0) + 

> gainma0*SPEED_0F_LIGHT*(l-l/2/gainma0/gainma0) 

*sin (phi) *Vector3 (0,1,0) ; 

> 

D.7 Copyright 

Copyright (C) 1994-2002 The Regents of the University of Cahfornia (Regents). 
All Rights Reserved. 

The code XOOPIC is referred to herein as the Software. 

Permission to use, copy, modify, and distribute the Software and its documen- 
tation for educational and research purposes, without fee and without a signed 
licensing agreement, is hereby granted, provided: (1) that the above copyright no- 
tice, this paragraph and the following two paragraphs appear in all copies, modifi- 
cations, and distributions; (2) you will not charge more than the cost of duplication 
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for copies of the original or derivative versions of the Software; (3) any export of 
the Software must be in comphance with U. S. export control regulations. For a 
license permitting for-profit distribution of the Software or its derivatives, contact 
The Office of Technology Licensing, UC Berkeley, 2150 Shattuck Avenue, Suite 
510, Berkeley, CA 94720-1620, (510) 643-7201. 

IN NO EVENT SHALL REGENTS BE LIABLE TO ANY PARTY FOR DI- 
RECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAM- 
AGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS 
SOFTWARE AND ITS DOCUMENTATION, EVEN IF REGENTS HAS BEEN 
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 

REGENTS SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUD- 
ING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANT- 
ABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE 
AND ACCOMPANYING DOCUMENTATION, IF ANY, PROVIDED HEREUN- 
DER IS PROVIDED "AS IS". REGENTS HAS NO OBLIGATION TO PRO- 
VIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MOD- 
IFICATIONS. 

For custom modification and support of the Software, contact Prof. J. P. 
Verboncoeur (johnv@eecs.berkeley.edu). If you use the Software for publication 
or other form of communication, as a courtesy, please use the following or similar 
citation: 

J. P. Verboncoeur, A. B. Langdon and N. T. Gladd, Comp. Phys. Comm. 87, 



199 (1995). Code available via |http://ptsg . eecs . berkeley. edu 



Appendix E 
OOPIC Input Files 

E.l A Sample Input File 

In this section a valid OOPIC input file is shown which can be used to recover the 
results presented in chapter El The values for (, vth and 7 can be changed easily. 

E-Layer 
{ 

Simulation of E-Layer 
} 



Variables 
{ 

c = 2.99792e8 
re = 2.82e-15 
e = -1.6e-19 
me = 9.11e-31 
muO = 1.256e-6 
epsO = 8.8542e-12 
gamma = 30 
zeta = 0.020 
vth = 0.002 
rO = 10.0 
L = 10*rO 
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range =20.0 

sigma = vth*rO/sqrt (l-2*zeta) 
beta = 1 - l/(2*gamma''2) 

Bext = 9. Ile-31*c*beta*gainma/1.609e-19/r0*(l+2*zeta) 
NO = gaiifflia*L*zeta/re 
dN = 0.0 
NMacPart = 5000 

norm = sqrt(2*PI)*sigma*2*PI*rO 
KK = NO/norm 
Q = NO*e 

I = Q*c/(2*PI*rO) 

} 

Region 
{ 

Grid 
{ 

J = 512 

xls = -range 

xlf = range 

K = 512 
x2s = -range 
x2f = range 
Geometry = 1 
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} 

Control 
{ 

dt = 5.0e-12 

BSinit = muO*I/L/2/PI*( 

EllipticE(2*sqrt (sqrt (xl~2+x2~2) *rO) * (sqrt (xl"2+x2~2) +rO) / 

( (sqrt (xl"2+x2''2)+r0) "2+r0''2*vth''2) ) * (sqrt (xl''2+x2''2) -rO) / 

( (sqrt (xl"2+x2''2) -rO) "2+r0~2*vth~2) - 
EllipticK(2*sqrt (sqrt (xl~2+x2~2) *rO) * (sqrt (xl~2+x2~2) +rO) / 

( (sqrt (xl"2+x2''2)+r0) "2+r0''2*vth''2) ) / 

(sqrt (xl"2+x2"2)+r0) ) 
BOS = Bext/L 
ElectrostaticFlag=0 

} 

Species 
{ 

name = electron 
m = me/L"2 
q = e/L 

collisionModel = 1 

} 

Load 
{ 

speciesName = electron 
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LoadMethodFlag = 

analyticF = KK*exp(-(xl''2+x2'^2-2*r0*sqrt(xl''2+x2''2)+r0"2)/ 

2 . 0/sigina'^2) * ( l+dN*sin(50*ataii2 (xl , x2) ) ) 

np2c = NO/NMacPart 

vldrift = c*(l-l/(2*gaiiima''2)) 

v2drift =0.0 

vSdrift =0.0 

xlMinMKS = -range 

xlMaxMKS = range 

x2MinMKS = -range 

x2MaxMKS = range 

} 
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